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ABSTRACT 






This report presents in detail a calculation method to compute the 
trajectory and ablation characteristics at the stagnation point of a 
glass sphere entering into the atmosphere of the earth from an arbitrary 
point in space. The underlying equations employed by the method include 
the transient effects, internal radiation, melting and nonequilibrium 
vaporization of the glass-liquid layer. A computer program written in 
Fortran IV language and a detailed description of the preparation of 
input for this program are included. The program is particularly applic 
able to the study of tektites and their atmospheric entry. The method 
and program could easily be altered to calculate the ablation at the 
stagnation point of a spherical glass tip on a missile-shaped body. 
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DEFINITION OF SYMBOLS 

(The following table lists the symbols that appear in the equations and 
the corresponding symbol of the computer output. The system of units 
is the meter-kilogram-second system where kg stands for kilogram mass . ) 


Equation 

Symbol 

Program 

Output 

Symbol 

Dimens ions 

Description 

Ai,A 2 ,A 3 

a 4 

t 

A1,A2,A3 

A4 


constants in vapor pressure 
function 

a 

VTTl 


m 2 sec/kg 

evaporation resistance 

a oo 


m/sec 

speed of sound 

b 4 

Bl,B2,B3 

B4 


constants in viscosity function 

C w,eq. 



equilibrium mass fraction of the 
injected vapor at the wall 

°P 

CP 

kcal/kg°K 

specific heat at constant pressure 
of the body material 

C D 

CD 


drag coefficient, Cj) = D/(p 00 W 2 A/2) 
where A is a reference area 

A = icr? 

r o 

^CO 

GINF 

m/sec 2 

gravity constant (free stream) 

h e 

HE 

kcal/kg 

enthalpy of air at the outer edge 
of boundary layer 


HW 

kcal/kg 

enthalpy of air at the wall 

*V 

HV 

kcal/kg 

heat of vaporization of the body 
material 

H 

H 

m 

geometric flight altitude 

% 


m 

flight altitude where slip flow 
regime begins 


v 







DEFINITION OF SYMBOLS (Continued) 


Equation 

Symbol 

Program 

Output 

Symbol 

Dimens ions 

Description 

k 

K 

kcal/m°K sec 

thermal conductivity of the body 
mater ial 

k air 


kcal/m°K sec 

thermal conductivity of air 

Km 

KM 


non-dimensional velocity gradient 
at outer edge of boundary layer 
dU e 

„ . 3T <*«„> 

“Sn W 

m 

MASS 

kg 

mass of the body 

mi 

Ml 


constant in vapor pressure 
function 

m o 

MO 

kg 

initial mass of the body 

M 

M 

! 

kg /kg mole 

molecular weight of air (free 
stream) 

^vap 

MVAP 

kg/kg mole 

; 

molecular weight of the gas 
vaporized 

M 

00 

MACH 

1 

! 

flight mach number 

n 

N 


the refractive index 

Nu 



Nusselt number 

P e 

PE 

kg/m sec 2 

pressure of air at the outer 
edge of boundary layer 

Poo 

PINF 

kg/m sec 2 

pressure of air (free stream) 

p 

vap 

PVAP 

kg/m sec 2 

vapor pressure of vaporizing gas 







DEFINITION OF SYMBOLS (Continued) 


Equation 

Symbol 

Program 

Output 

Symbol 

Dimens ions 

Description 

P* 

vap 

PVAPS 

kg/m sec 2 

equilibrium vapor pressure of 
vaporizing gas 

q 

aero 

QARO 

kcal/m 2 sec 

aerodynamic heating rate for zero 
vaporization 

q rad 

QRAD 

kcal/m 2 sec 

radiative heat flux rate from the 
body wall 

o 

u 

RFO 

m 

initial radius of the body 

r f (t) 

RFT 

m 

instantaneous radius of the 
spherical body cap 

R e 

RE 


Reynolds number based on initial 
body diameter 

^ef f 

R(EFF) 


effective reflectivity of the 
surface 

t 

TIME 

sec 

time 

At 

DT 

sec 

time grid in difference method 

T 

T 

°K 

temperature 

St/Sy 

TP 

°K/m 

temperature gradient in Y-direc- 
tion 

S 2 T/Sy 2 

TP 2 

CM 

0 

second derivative of temperature 

St/ St 

DTDT 

°K/ sec 

derivative of temperature with 
respect to time 

T, 

oo 

TINF 

°K 

free stream temperature 

Tw 

TW 

°K 

surface temperature of body 


vii 








DEFINITION OF SYMBOLS (Continued) 


Equation 

Symbol 

Program 

Output 

Symbol 

Dimensions 

Description 

T o 

TO 

°K 

temperature of body at initial 
time 

T e 

TE. ' 

°K 

temperature of air at outer edge 
of boundary layer 

u 

UINF 

m/ sec 

horizontal free stream velocity 
component 

dU/dt 

DU/Dt 

m/sec 2 

derivative of flight velocity 
component with respect to time 

V 

VINF 

m/sec 

vertical free stream velocity 
component 

dV/dt 

DV/Dt 

m/sec 2 

derivative of flight velocity 
component with respect to time 

V 

V 

m/sec 

ablation rate 

v w 


m/ sec 

vaporization rate v w - v(Y=0) 

V 


m/ sec 

total ablation rate at each time 

00 



Voo = v(Y-Y 0 ) 

w 

WINF 

m/sec 

free stream flight speed 

dW/dt 

DW/Dt 

m/sec 2 

acceleration 

w c 


m/sec 

circular speed of the earth 

Y 

Y 

m 

coordinate measured from body 
surface 

£Y 

DY 

m 

distance grid in difference 
method 


viii 







DEFINITION OF SYMBOLS (Continued) 


Equation 

Symbol 

Program 

Output 

Symbol 

Dimens ions 

Description 

y b 

YB 

m 

' upper limit of integral in heat 
balance equation 

Yo 


m 

point in the body where no melt- 
ing occurs 

Y s 

YS 

1 m 

1 thickness of body material lost 
due to ablation 

Y s 

w 

YSW 

m 

thickness of body material lost 
due to vaporization 

a 

ALF 

1/m 

reciprocal radiation mean free 
path 

°A 

ALFA 

l/m 

absorption coefficient 


ALFV 


vaporization coefficient 

e 

E 


emissivity constant of the 
opaque body material 


PHI 

degrees 

angle of attack measured from 
horizontal 

M 

MU 

kg/m sec 

viscosity of the body material 

Me 

1 

MUE 

kg/m sec 

viscosity of air at outer edge 
of boundary layer 

^00 

MUINF 

kg/m sec 

viscosity of air (free stream) 

p 

RHO 

kg/m 3 

density of the body material 

Poo 

RINF 

kg/m 3 

density of air (free stream) 

Pe 

RHOE 

kg/m 3 

density of air at the outer 
edge of boundary layer 






DEFINITION OF SYMBOLS (Continued) 


Equation 

Symbol 

Program 

Output 

Symbol 

Dimens ions 

Description 

V* 

TAW 

kg/m 2 sec 3 

shearing stress at wall divided 
by X 

t 

SI 


heat blockage factor, a cor- 
relation function derived from 
solutions to the boundary layer 
equations 


INT 

°K m/sec 

value of integral in heat 
balance equation 

Px 

BETAl 


constant in shear stress relation 

Sf/Sy 

DF 


see equation (A-ll) 
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A CALCULATION METHOD FOR THE ABLATION OF GLASS-TIPPED BLUNT BODIES 

SUMMARY 


This report presents in detail a calculation method to compute the 
trajectory and ablation characteristics at the stagnation point of a 
glass sphere entering into the atmosphere of the earth from an arbitrary 
point in space. The underlying equations employed by the method include 
the transient effects, internal radiation, melting and nonequilibrium 
vaporization of the glass-liquid layer. A computer program written in 
Fortran IV language and a detailed description of the preparation of 
input for this program are included. The program is particularly applic- 
able to the study of tektites and their atmospheric entry. The method 
and program could easily be altered to calculate the ablation at the 
stagnation point of a spherical glass tip on a missile-shaped body. 


I. INTRODUCTION 


A great many investigations have been conducted in the past eight 
years to determine the ablation of bodies entering the earth's atmosphere 
at hypersonic speeds. For glassy or quartz-like material, the ablation 
follows a process of both melting and vaporization. Sutton [15] in 1958 
obtained the first exact numerical solution for the steady state abla- 
tion of glass-like material at the stagnation point of a body subjected 
to hypersonic flight conditions. His method is limited, however, in 
that it accounts for melting only, and even for this restricted case, it 
is very difficult to employ. Bethe and Adams [2], Few and Fanucci [10], 
Scala [14], and many others presented quasi-steady solutions that accounted 
for both melting and vaporization. Adams [1] presented a calculation 
method that employed a finite difference procedure. This procedure 
accounted for the transient effects of the ablation problem which 
included both melting and vaporization effects. Scala and Vidale [16], 
who considered the nonequilibrium effects of vaporization in the quasi- 
steady approximation found that, when a material is subjected to severe 
heating conditions, the net rate of ablation due to vaporization may be 
diffusion controlled, kinematically limited, or both. Kadanoff [9] 
developed a method for calculating the net radiative flux within a semi- 
infinite body which emits, absorbs, and scatters this radiation and which 
allows some radiation to escape from its surface. Chapman [3] in 1963 
presented a condensed account of a nonsteady calculation method which 



accounts for melting, vaporization, and internal transport of radiant 
energy. Chapman includes approximately the effects of the presence of 
O 2 and hence production of the S^O vapor as given by Hidalgo [7]. 

A finite-difference method is presented for the solution of the 
ablation problem in the vicinity of the stagnation point which accounts 
for transient effects, nonequilibrium vaporization, and internal radia- 
tion for a spherically or hemispherically shaped body composed of a 
glassy or quartz-like material that melts and vaporizes. The equations 
developed by Chapman [4] for entry into a planetary atmosphere are 
employed in this calculation scheme to compute the trajectory of the 
body entering the earth 1 s atmosphere. The equations given in the follow- 
ing sections of this paper are a system that must be calculated for each 
time step and are presented in the order of computation in the scheme. 

At each time step, there is an iteration involving T w (the surface 
temperature) such that a heat balance equation at the surface is satis- 
fied. 


A computer program written in Fortran IV language for the equations 
in this report is presented in the appendix. This program is designed 
to calculate the trajectory, temperatures, and ablation history at the 
stagnation point of a sphere or hemisphere composed of a glassy material 
which is entering into the atmosphere of the earth. The preparation of 
16 input data cards is required for each particular example to be com- 
puted. The cost of running this program on the 7094 computer at Marshall 
Space Flight Center is approximately one cent per time step for the case 
without internal radiation and 2.5 cents per time step for the case with 
internal radiation. 

The author is indebted to Mr. Verkuel Eubanks of the Computation 
Laboratory, General Electric Company, Huntsville, Alabama for the detailed 
programming and valuable suggestions concerning the numerical solution of 
this problem. 


II. INITIAL CONDITIONS AND CALCULATION PROCEDURE 


The entry into the earth's atmosphere begins at an initial altitude 
Hq (an input value, usually 150,000 m) ; i.e., t = o when H = H 0 . Other 
initial conditions which must be given are angle of attack, 0 O , initial 
flight velocity, W Q , radius of the sphere or hemisphere, rf Q , and the 
physical properties of the body. The preparation of these input data 
for the computer program is outlined in Appendix B. At the initial time 
it is assumed that the body has a uniform temperature T 0 = 300°K through- 
out the body. 



t 


The differential equations describing the flow of the viscous glass 
layer in the vicinity of the stagnation point are given in Appendix A 
of this paper. The mathematical procedure solves each of these three 
basic conservation equations by taking small step-by-step increments in 
time. The equations presented in the following sections are a system 
that must be calculated for each time step; it is assumed that the cal- 
culations are being carried out at time t and that the solution is known 
at time t - At. In order for convergence of the forward difference pro- 
cedure in the solution of the transient energy equation, the following 
relationship between grid and material properties must be satisfied: 


A 




where the material properties p, Cp, and k are assumed constant and not 
a function of temperature. If the above relationship is not satisfied, 
the program will solve the equation with the equality sign for AY, and 
the calculation will continue. 

The computer program contains an option of whether to account for 
or disregard the effects of internal radiation. The first page of the 
computer output is printed for the appropriate case as: 

(a) THE ABLATION PROGRAM WITH INTERNAL RADIATION 

(b) THE ABLATION PROGRAM WITHOUT INTERNAL RADIATION. 

The case without internal radiation assumes the body to be composed of 
an opaque material, and only heat radiated away from the surface is 
accounted for. The case with internal radiation uses the equations 
derived by Kadanoff [9]. 

Calculation by the computer program ends when one of the following 
conditions are met: 

(1) The body reaches the surface of the earth (H ^ o). 

(2) The body is dissipated due to ablation. 

(3) The altitude H exceeds a predetermined maximum. This is 

to avoid a body from bouncing out of the earth's atmosphere 
and continuing indefinitely. 

(4) The flight time reaches a predetermined limit. 

(5) An error condition exists. 


3 



In each case an appropriate message is printed explaining the reason for 
terminating the calculation. 
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The program allows for up to thr^e different time steps and print 
frequencies during any one run. For instance, let At be the delta time, 
T m be the upper time limit at that time step, and Mp be an integer 
denoting the number of time steps between prints. 


Time Step 


Maximum Time Print Frequency 


At x 


T 

m x 


M 

Pi 


At 2 



M„ 


2 


At 3 



M 


P 3 


The program would use a At of At x between time = 0 and time = Tm^ print- 
ing every Mp^ steps. The program would then use At 2 between T mi and Tm 2 > 
etc. The program selects the largest of the Tm values as the maximum 
time for the entire run and selects the largest At values to check the 
grid ratio discussed above. If one wants to hold the print interval and 
At constant throughout a run, it is best to set Tm x equal to end of 
flight time and Tm 2 and Tm 3 equal to 0. 

The subsequent sections, which give the equations used in the cal- 
culation method, are presented in the order that the computer program 
actually solves the problem. Since most of the underlying analysis of 
the equations is well known, they are presented in their final form. 

Many of the relations were taken from other references dealing with 
aerodynamics, ablation, and boundary layer theory. 

The computed results by this program for a typical atmospheric 
entry of an initially spherical-shaped tektite are shown in Figure 1. 

The initial flight conditions and physical properties for this example 
are 


(1) case with internal radiation 

(2) H 0 = 150,000 m 

(3) W Q = 11,200 m/sec 


(4) 0 O = -20° 


4 


t 


(5) r £ = .01 m 

o 

(6) k = 3.75 x 10“ 4 kcal/ (m °K sec) 

(7) p = 2400 kg/m 3 

(8) c p = .34 kcal/ (kg °K) 

(9) by = 3050 kcal /kg 

( 10 ) \ ap = 40.278 kg/kg-mole 

(11) constants in viscosity function (see equation (78)) 

Bi = .1 
B 2 = 27,620 
B 3 = 262 
B 4 = -9.09 

(12) constants in equilibrium vapor pressure function (see 
equation (80)) 

Ai = 101,325 

A 2 = 0 

A 3 = -57,800 
A 4 = 19.1 

(13) = .28 

(14) 0 ^ = 0 o 

(15) a = 1900 1/m 

(16) a A = 285 1/m 

(17) n = 1.5 
(IS) Reff = »2 


5 









III. RADIUS OF CURVATURE OF THE FRONT FACE OF THE BODY 


At the initial time, the body is assumed to be either a sphere or 
hemisphere. During the entry of the body into the atmosphere, the 
radius of curvature rf(t) varies with the extent of ablation Y s (t). 

For bodies where the ablated thickness Y s (t) is small compared to the 
initial radius of curvature ff(t) such as the ablation of the front 
face of a blunt missile body, the variation of r-p(t) can be neglected 
in the calculation scheme. For small bodies such as tektites where 
the ablation thickness is comparable to the body radius, the variation 
of T£ (t) is important and therefore must be determined at each time 
step. Two methods of determining this radius of curvature are given 
below. One method (Method 1) gives results which overestimate the 
total mass lost and the other method (Method 2) probably underestimates 
the total mass lost. The program is coded such that either method can 
be used in the calculation scheme. 

A. Method 1 

This method assumes that the mass lost by the body is the mass 
lost due to both melting and vaporization; i.e., the material that is 
melted and vaporized is removed from the body. The assumed body shape 
is shown in the following illustration (Figure 2). 



Figure 2. The Body Shape for Method 1 
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The radius of curvature is given by 


r f (t) 


- r , 


o 



Y 


2 

S 



(i) 


where the total thickness along the axis of symmetry that is lost up 
to the time t due to the ablation is 


Y s (t) - - 


dt 


o 


( 2 ) 


and VooCt) is the instantaneous ablation velocity. When the condition 
Y s (t) ^ Vf Q is met, the computer program is coded such that the calcula- 
tion stops and a message is printed out denoting that half (all) of the 
sphere (hemisphere) has been dissipated due to ablation. 

B. Method 2 


Chapman [3] in his experimental studies on the ablation of a 
tektite glass sphere placed in a hypervelocity arc jet, found that the 
melting begins at the stagnation point where the heating is most severe 
and the molten material flowed around the sphere and solidified on top 
of the original spherical surfaces. The mathematical model for this 
body shape is shown in Figure 3. It is assumed that the mass loss of 
the body is the mass lost due to vaporization only and that the mass 
melted is equal to the mass of the flanges. Chapman derived the follow- 
ing empirical relation for the radius of curvature rf(t) as a function 
of ablated thickness Y s (t) from the experimental results: 



o o 


The thickness of body material that has been lost due to vaporization 
is 



(4) 
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where v w (t) is the instantaneous vaporization rate at the surface of 
the body at the stagnation point. 



C. Calculation of the Body Thickness 

The thickness of the body S(t) measured along the axis of 
symmetry is given at each time by 

S(t) = L(t) • AT. 


9 




The initial thickness of the body is given by S Q = L Q • AY. 


L 


o 


i • rf 
L o 

AY 


where i = 1 for a hemisphere and i = 2 for a sphere. The computer 
program approximates the function L(t) by taking it as the integer part 
of the number 


S - Y (t) 
o s 

AT 


+ 


.5. 


Therefore, L(t) + 1 is the number of points Y along the axis of symmetry 
that are considered in the calculation at time t. 


IV. THE TRAJECTORY OF THE ENTRY BODY 


By neglecting the lateral forces. Chapman [4] studied the descent 
of a body in a meridian plane of a spherically symmetric atmosphere 
about a spherically symmetric planet. The coordinate system and 
velocity components are shown in Figure 4. 
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Figure 4. Coordinate System for Trajectory Calculations 






For a spherical body the resulting equations of motion are as follows: 


dU(t) 

dt 


d . V . 0 0 

dt 


where 


-U(t) V(t) 


H(t) + R . . m(t) 

v ' sea level 


P (t) 

C D (t) ( U 2 (t) + V 2 (t) 


+ .75 pv (t) 


U(t) JtTj 


JvHt) + V 2 (t) 


(5a) 


= " 8 m (t) + 


u 2 (t) 


rpJ fc ) 


«> w H(t) + R . . 

' ' sea level 


2 D' 


C (t) ( U 2 (t) + V (t) 


+ .75 pv^(t) 


V(t) ki ‘ 


m(t) s/u 2 (t) + V 2 (t) 


(5b) 


(a) H(t) = H(t - At) + V(t) . At is the local geometric altitude. 


g o 0 ( t ^ ®sea level 


+ R sea level N 2 . 


R 


sea level 


is the local value 


of gravitational acceleration. R se a level - 6,380,000 m and 

e = 9.80665 m/sec 2 for the system of units used herein, 

sea level 

(c) p (t) = the density of air in free stream, a function of H(t). 

(d) p = the density of the body material 

(e) Cjj(t) = the drag coefficient. 


11 



The ordinary differential equations (5a) and (5b) are numerically solved 
by the method of Runge-Kutta at each time t for the flight velocity com- 
ponents U (t) and V(t). To avoid including the above differential equa- 
tions in the iteration procedure for the surface temperature T w (t), it 
is assumed that, in (5a) and (5b), v w (t) = v w (t - At) and m(t) = m(t - At) 
which is a good approximation due to the smallness of the time step At. 


Relations for the drag coefficient Cp(t) for a sphere and a hemi- 
sphere were derived by empirically fitting curves from various sources. 
These formulas for different flight conditions and regimes (see Sec- 
tion VIB) are as follows: 

(a) Molecular Flow Regime (H > ftp) for Spheres and Hemispheres 

(1) 0 § M ^ 9; C_ = + 1.68 

v 7 oo ’ D M 

00 

(2) 9 <M oo ^ oo; C D = 2.0. 

(b) Continuum Flow Regime (H < ftp) for Spheres 

(1) 0 S M ro S .8; C D = .5 

(2) .8 <M co S 1.26; C D = .812 - .023 

(3) 1.26 < ^ 2.0; C D = 1.034 - .027 

(4) 2.0 < i co; C D = .9 + M 0O /'/Ri • 


(c) Continuum Flow Regime (H < Hp) for Hemispheres 


(1) 0 S Mco S 2.0; C D 

(2) 2 < S oo; C D = 


= 1.35 

1.35+ M^A/R^ . 
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The free stream Mach number and Reynolds number are given by 


M (t) 


:= W(t) 


a (t) 

00 


( 6 ) 


Re< t ) = 


p ro (t) W(t) 2r fo 


(7) 


where the atmosphere properties p^, a^, T^, and are given as a func- 
tion of the flight altitude, H. The computer program employs a sub- 
routine that contains a table of atmospheric properties taken from 
Reference 17. This subroutine interpolates the table for the properties 
at the given value of H for each time. After solution of the two dif- 
ferential equations (5a) and (5b), the following relations are calculated: 

(a) The resultant velocity 


W(t) = Vu 2 ( t) + V 2 (t) . 


( 8 ) 


(b) The angle of attack (see Figure 4) 


0(t) = tan -1 


V (t)/U(t) . 


(9) 


(c) The acceleration of the body 


Sw 

dt 


•£ + 




w 


( 10 ) 
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V. PROPERTIES OF AIR BEHIND THE NORMAL SHOCK 


Stagnation point properties of air behind the normal shock, or at 
the outer edge of the air boundary layer, that are needed in the abla- 
tion problem are approximated by curve fits of curves and tabulated data 
for air in thermal and chemical equilibrium published in Reference 8. 

The properties and their respective curve fits are as follows: 

A. Temperature 

(1) For W > 2100 (m/sec) and M^ < 35, 


T e , 


ideal 


T°o ideal 


where 


T e , 


M c 


ideal , oo 

T 5 * 


The ratio of real to ideal gas effects is given by 


(ID 


( 12 ) 


d Q + d^ + d^ + d^ + d^£ + dsMcL 

e ideal 


(13) 


where 

d Q = 4.016949 - 1.49287 x 10" 5 H 
di = -.895475 + 4.51127 x 10 -6 H 
d 2 = 9.28796 x 10" 2 - .54521 x 10 -6 H 
d 3 = -4.746323 x 10 -3 + 3.05088 x 10 -8 H 
d 4 = 11.6111955 x 10“ 5 - 7.96241 x 10 _1 ° H 
d 5 = -10.86105 x 10" 7 + 7.84415 x 10" 12 H 

and H = H(m) is the flight altitude. 
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(2) For W § 2100 (m/sec) and M ra < 35 


T e 

Te ideal 

(3) For ^ 1 35, 



B. Density 

(1) For W > 2100 (m/sec) and < 35, 
p e p e Pe ideal 


where 


P °° Pe ideal 


'ideal 


is given by one of the two following relations depending on Moq: 


(a) 


e ideal 


6M 2 - 7 / 2 

00 

M 2 + 5. 


7M 2 - 1 


W 


1 + 


5/2 


for Mo,, - 1* 


(b) 


for M < 1. 


( 14 ) 


(15) 


(16) 


(17) 


(18) 
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The relation for the ratio of real to the ideal gas value is given by 


Pe ideal 


" a c + a i Moo + a 2 


M 2 + a 3 

00 


K + 


I 4 M 


00 


where 


a 0 = .269623 - 2.47746 x 10" 2 S 
a x = .1975914 + 7.249941 x 10“ 3 S 
a 2 = -1.334415 x 10* 2 - 6.96885 x 10“ 4 S 
I 3 = 5.06022 x 10" 4 + 3.2325 x 10" 5 S 
a 4 = -.832101 x 10* 5 - 6.34127 x 10" 7 S 
a 5 = 3.569058 x 10~ s + 3.77445 x 10" 9 S 

and S = H/1000. 

(2) For W ^ 2100 (m/sec) and M m < 35, 

_^_ =1 

Pe ideal 


(3) For ^ S 35 



C. Pressure 


(1) For W > 2100 (m/sec). 


P P e . 

o c “ 


ideal 


P Pp, 


00 e ideal oo 


( 22 ) 


where 


Pe, 


ideal 


is given by one of the two following relations depending on Mo,: 




(b) 


e ideal 


M c 


7/2 


1+-) for < 1, 


(24) 


The relation for the ratio of real to the ideal gas value is given by 


Pe, 


ideal 


C Q + C ± M* + C 2 + C 3 + C 4 + C 5 Jfg, 


(25) 


where 


C c = 1.0099865 - 2.722132 x 10" 6 H 
C x = -4.9934 x 10" 3 + .83986 x 10 -6 H 
C 2 = 19.39623 x 10 -4 - .860322 x 10 -7 H 
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VI. THE SOLUTION OF THE ABLATION PROBLEM 


A. The Iteration Procedure to Determine the Wall Temperature 

The surface temperature, T w , is determined by an iteration 
procedure which is satisfied when the heat balance equation at the sur- 
face-air interface is satisfied. The first guess at time t for T w is 


T w 1)(t) = v* - + • At * (31) 


It is desired to find a value of T w which will satisfy the heat balance 
equation (67). After at least two values of T w are tried and do not 
satisfy the relationship, one can get a good approximation of the value 
of T w by plotting T w vs E, where E is the error in equation (67), and 
finding the point where the error would be zero. 



It is not necessary that the error be of opposite sign, not is it neces- 
sary to use a higher order interpolation, since the first guess will be 
very good and the second guess will yield an E of opposite sign or closer 
to zero. This can be done because the left-hand side of (67) monotonically 
Increases and the right-hand side of (67) monotonically decreases with an 
increase in the size of T w . By putting a straight line Tw - a + bE which 
goes through the points E x and E 2 , we can get the next guess for T w by 
evaluating the equation T w = a 4 - bE for E = 0. This reduces the problem 
to finding "a” since T w = a for E = 0. Using Kramer’s rule, 


T Wl ~ a + bEj 
Tw 2 = a + bE 2 


a 
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For subsequent iterations, the last two points will be used. Care must 
be taken to avoid letting E 2 - Ei become too small, causing overflow. 

The program includes a test such that a guessed value of T w would never 
yield a vapor pressure Pvap(Tw>t) larger than the pressure in the bound- 
ary layer P e (t). 

B. Calculation of the Aerodynamic Heat Flux 

For given values of the surface temperature, flight altitude, 
and flight velocity, relations are given for the aerodynamic heat flux 
to a nonvaporizing wall q aero for the different flight regimes. The 
flight regimes, in the order in which an object entering the earth’s 
atmosphere encounter, are the free molecule, transition, slip, and con- 
tinuum. It was found that formula (32), given below for the aerodynamic 
heat flux employed by Reference 3 was a good approximation from the 
initial time when the object is in the free molecular regime through the 
transitional regime. The altitude, Hj, at which the object enters the 
slip flow regime, i.e., passes from the transition regime to the slip 
flow regime, is assumed to occur when the Knudsen Number (Kn) reaches a 
value of 0.1. Figure 5 presents the altitude Ht as a function of the 
body radius for the earth's atmosphere. 
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Figure 5. Altitude Where Slip Flow Regime Starts 
As A Function of Nose Radius 




The aerodynamic heat flux q erQ for the slip flow regime is 
approximated by a smooth transition from the transitional regime, equa- 
tion (32), to the continuum flow regime, equation (34). This transition 
takes place over a prescribed number of time steps beginning at the time 
when the altitude H =^11^. 

Relations for the aerodynamic heat flux taken from various 
sources for the different flight regimes are as follows: 

j(l) For 150,000 ^ H(m) ^ Kp - Reference 3 found the following 
equation to be a good analytical representation of the experimental data 
of Reference 6 in the molecular and transitional flow regimes: 


q aero ^aero 

- m °l cont _ - 

q a«° ’ As 5 Tp — ” V ro ’ 

/ H aero aero Tran 

\ / mol cont 


(32) 


where q aero is the aerodynamic heat flux for the free molecule regime 
mol 

given in Reference (11) for a sphere as 


q 

n aero 

mol 


.53733 x 10“ 4 


p W 

^00 



1.9262 x 10+ 3 (T w - T J 


(kcal/m 2 sec). 


(33) 


where p = p (kg/m 3 ), W = W (m/sec), and T = T (°K). q aero is the 
00 00 cont 

aerodynamic heat flux in the continuum flow regime given by Reference 5 
as 


q 

n aero 

cont 


23,812.9 


/ 


r f (t) 


(w/w c ) 3 * 15 (r 

\ e 



(kcal/m 2 sec). 


(34) 
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where 


h = 

e 8374.88 


+ .24 T 


(kcal/kg) 


(35) 


I 397.92346 x 10 xa 
■J 6.37 x 10 B + H 


(m/sec) , 


(36) 


and the enthalpy at the wall for air is approximated by 


h = -5.3944 + .24019 T w + 2.03371 x 10* 5 T 2 - 2.30515 x 1(T 9 T 3 

W W w 

+ .96954 x lO" 13 (kcal/kg). (37) 


(2) Transition Regime from Slip Flow to Continuum Flow - At 
the first time t, when H ^ H T which is designated t*, the following 
transition equation for q aero is introduced. It is completed in ten 
time steps; thus, ten time steps after H ^ By, it is assumed that the 
re-entering object is in continuum flow. Let t* + 10 At = t*. 



" 


1 

* H 

U 

1 

q - q + 

n aero aero 

4 

aero 

4 

aero 


ft 

J * 

1 

rt 

* 

tran 

cont 

tran 


L 1 


(38) 


(3) From the time tg until time t^ (t’3 is defined as the 
first time that the relation W = 2100 (m/sec) is violated), the con- 
tinuum flow relation for is used; i.e., 

d. ci O 


q 


aero 


q 

aero 

cont 
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(4) From the time t* until impact time, the following rela- 
tion for the aerodynamic heat flux which was taken from Reference 13 
is used: 


(.715 )- 4 yip^ <3,) 

(kcal/m 2 sec). 

The thermal conductivity of air at the wall 1^ a £ r i- s approximated by 

k . = .6716646 x 10" 8 + .2429834 x 10" 7 TL. - 1.811997 x 10” 11 T 2 

w,air ** w 

+ 1.3873689 x 10" 14 T 8 - .5989437 x 10" 17 T^f + 1.0229805 x 10" 21 

[kcal/(m°K sec)], (40) 


aero 




(T e 


r w) 

V w 


and the viscosity of air at the wall a ^ r is given by the Sutherland 
law: 

VL . 14 ‘ 7 ?in 3 x 10 ” 7 ^ kg ^ m sec ^* ( 41 ) 

'w.air v w . , 110 

Tw 

The relation for (NuA/Ri)^^ is approximated by the following curve 
fits of results in Reference 13; 

T 

(a) For < 1.6, 

L e 

T 

(Nu/VRi) pr=1 = .705 + .055 (42) 
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(b) For ~ ^ 1.6, 
1 e 


(Nu/N/Re ) pr=1 


.755 + .025 . 

T e 


( 43 ) 


C. Heat Blockage Factor 

Equations (44) - (48) in this and the next section are applic- 
able to the continuum gas dynamic regime. In the free molecule flow 
regime, there is no heat blocking effect by the vaporizing species; 
i.e., \|/ = 1. These equations were employed at all times in the com- 
puter program; however, no significant vaporization will generally 
result before the continuum flow regime is entered. The heat blockage 
factor, according to Reference 2, due to vaporization is 


* = 


1 


1 


C 

w,eq. 



.68 (M/M^p)* 26 


(44) 


where C w , the equilibrium mass fraction of the injected vapor at 

the wall’ is given by 


C 

w,eq. 


1 


1 + 


M 


M 

vap 




(45) 


M is the molecular weight of air, a function of altitude, and M va p is 
the molecular weight of the vaporizing gas. The equation for the vapor 
pressure P V ap is given in the section on material properties. 


D. Ablation Rate at the Wall 


The ablation rate at the wall, v w , is due to the vaporization 
process only because the melting component of ablation is zero at the 
wall. By combining the boundary layer solution for given in Refer- 
ence 2 for a Lewis number of one with Scala f s [16] kinetic theory 
expression for C^, the following equation for v w results: 
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where a-vm is the resistance of the material to the vaporization process. 
It is given by Scala [16] as 


J 

T 

vap w 

a 

V 

(P-P )M + P M 

e vap vap vap 


(47) 


where (universal gas constant) = 8.314 x 10 3 kg m 2 / (kg mole sec 2 °K) 
and is the vaporization coefficient. When a^ = 0, the following 
equation for v w results: 


v = 
w 



(48) 


The computer program interprets a zero input value for a v to mean 
a v = oo which due to (47) makes a-^ = 0; thus, when (^ == 0 is an input, 
v w is computed by equation (48). 


E. The Temperature Profile 

(a) The surface temperature, T^, is given by an iteration 
process. Section VI-A. 

(b) The forward difference procedures gives the temperature 
profile for 2z^Y § Y § (L - 1) AT: 


T(Y, t) = T(Y, t - Z£) + • Ztf, (49) 

where T.AY is the thickness of the body along the axis 
of the symmetric body. 
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The temperature at the last station L/^ is given by 


T(LZST,t) = T[(L - 1) • £f,t]. 


( 50 ) 


(c) The temperature at Y = is given by the following curve 
fit: 


T(£Y,t) = i T(o,t) + T(22Sf,t) - i T(32ST,t). (51) 


F. The Ablation Rate Profile 


After having a complete temperature profile, the ablation rate 
profile can be calculated by the following equation which results from 
the continuity, momentum, and the wall ablation rate equations (see 
Appendix A): 


Y Y, 


v (Y, t) = 


v (t) - 
w 


{; 


[ V‘> 


“ PocXt)] + 


2p 


r f (t) 


—I Y + 2 — 
dtj X 


p(Y, t) 


dYdY 


(m/sec), 


(52) 


where Y 0 is a point in the body where the integrand is approximately 
zero; i.e., the material is in a solid state at Y c . The viscosity of 
the material p(Y,t) is given in the section on material properties as 
a function of the temperature T(Y,t). 

The shearing stress relation is given by 


where 


~X = (T w /X) o W 1 + ^ )] > 


(P e /P va p) " 1 


(53) 


(54) 
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and p x = constant. The shearing stress at the wall for the nonvaporizing 
case, (%/X) 0 , is given for a sphere by Reference 12 as 


(t /X) = 

w ' o 


.727 


W 


L r f( t > 




5 


(kg/m 2 sec 2 ). 


(55) 


where p^(i = e,w) is given by the Sutherland law 


v -fh 


14.76303 x 10 


,-7 


1 + 


110 

T. 

l 


(kg/m sec). 


(56) 


G. The Energy Equation 

The derivative of the temperature with respect to time can 
now be calculated. 

(1) For Y = 0 and Y = 1 ST/St is approximated by a back 
wards difference quotient. 


ST(Y,t) = T(Y,t) - T(Y, t - At) 
St At 

(2) For 2/^Y SYS L£Y, 


(57) 


St(Y. t) = _k_ S 2 T(Y t) . ST(Y, t) 1 SF(Y.t) 

St pc p Sy^ (Y ’ ' Sy pc p Sy ’ 


(58) 


where the flux of radiative energy, F, is defined by equation (A-ll). 
For the case without internal radiation, SF/SY = 0. The second deriva- 
tives of T with respect to Y in equation (58) are given by the follow- 
ing relations: 
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(a) For Y = 0 and Y = 1AY, 


cl 2 T(Y,t) 


pc 


T(Y.t) - T(Y, t 
At 




( 59 ) 


(b) For 2ZY S Y ML - 1) ZY, 


d 2 T(Y, t) _ T(Y + £Y,t) - 2T(Y,t) + T(Y - AY,t) 

(AY) 2 


(60) 


and the first derivatives of T with respect to Y are given by the 
following: 


(a) For Y = 0, given in Section VI-H of this paper. 

(b) For 1AY ^ Y i (L - 1) AT, 


5T (Y, t) = T(Y 4- AT.t) - T(Y - AY, t) 
SY 2 AY 


(61) 


At the last grid point LAY, both the first and second derivatives of T 
with respect to Y are zero. 


H. The Heat Balance Equation 

A heat balance equation at the gas-liquid interface which 
must be satisfied at each time t is now derived, and is based on the 
fact that no heat can be stored at the gas-liquid interface. The net 
flux of heat on the gas side of the interface must equal the net flux 
of heat on the liquid side. Hypersonic entry into the earth's atmos- 
phere of a blunt-nosed object generates a curved detached shock wave 
which results in very large (gas-cap) temperatures. The radiation, 
termed here gas-cap radiation, from the high temperature air is approxi- 
mately proportional to the nose radius of the object. This heat addi- 
tion is neglected in the method presented herein since we were concerned 
only with objects of very small radius such as tektite and australite 
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bodies. However, gas-cap radiation should be accounted for when larger 
bodies are being considered. On the gas side, there are the following 
amounts of heat either arriving at or leaving the interface: 

(1) the aerodynamic heat flux (<i aero i|r) arrives at the inter- 
face, 

(2) the heat flux taken up by the vaporization process, 

( (q vap = "f^w M leavin S the interface, and 

(3) the heat flux radiated (q rad ) which leaves the interface. 
The heat fluxes on the liquid side of the interface are 

(1) the heat flux conducted at the interface, q c = -k(c)T/dY) w , 
and 

(2) the heat flux being radiated up to the interface from the 
interior of the body, (q rad ); this term is zero for the 
case without internal radiation. 

Equating the net heat flux on the gas side of the interface with the 
net heat flux on the liquid side yields 


^aero ^ %ap ^rad ^c ^rad’ 


(62) 


which can be written as 


q ilr + pv h = -k(£T/dY) . 

4 aero v H w v w 


(63) 


For the case without internal radiation, this equation is 


q \|f+pvh-q J = -k(c)T/c)Y) , 
^aero Y H w v rad w’ 


(64) 


where q rad is given by Stefan's law for emission of energy from the sur- 
face of a body: 


q , = e crT , 
H rad w 


(65) 


where a = 1.378 x 10” 11 [kcal/m 2 sec(°K) 4 ]. 
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The heat conduction term -k(ST/dY) w in equation (63) is deter- 
mined by integrating the energy equation (A-3) over Y from Y = 0 to Y = Y^ 

where Y B is any point within the body. It is best to make Yg at least 

equal to 6/ST to perform the numerical integration involved to an agree- 

able degree of accuracy. Integrating the energy equation (A-3) yields 



Combining equations (63) and (66) yields the following equation, which 
is subsequently called the heat balance equation: 


q xk + pv h = 
n aero Y w v 


- k 


dY 


+ pc 


clT . ciT , , 

at + v 1 dY + 


w 


(67) 


For the case without internal radiation the heat balance equation is 


q i 

n aero 


+ pv 


w 


h - 

v 


*rad 



( 68 ) 


The iteration for T w is ended when the heat balance equation is satis 
fied within a prescribed tolerance. 


I. Calculation of the Mass m(t) 

The equations presented in this section are based on the body 
shapes given in Figures 2 and 3. 

(1) Mass for Method 1 Geometry (see Figure 2) 

For this method, it is assumed that all of the material 
that melts and vaporizes is removed from the body. The equation for 
the mass of the body which varies with time due to the ablation proces- 
ses is 
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(69) 


m(t) = m o - P jt 


'! Y S 
L- O 


r y 2 Y 3 - 
r o s , _JL 

2 6 


where tn 0 is the initial mass of the body. 

(2) Mass for Method 2 Geometry (see Figure 3) 

For this method it is assumed that 

(a) mass lost = mass lost due to vaporization and 

(b) mass melted = mass of the two flanges. 

The points (| x t^) shown in Figure 3 are calculated by 


lx = 


r f D - r f < c > + hS 
2h 


and 


Til 


= + 


2 

r fn ' 


[rf o - r^(t) + h 2 ] 2 
4h^ 


where 


h - - r f„ + Y s 


and Y g is given by equation (2). The radius r^ is given by 


rf v = 2 (Sl + Sl) + 2(1! + 8l ) 


Sl rf o ' Y ®w' 


where Y„ is given by equation (4), 


J w 


(70) 

(71) 

(72) 

(73) 

(74) 
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( 75 ) 


j = r f - S-l. 
v 

The radius rf(t) is given by equation (3). The volume lost due to 
vaporization is given by 


AV vap - * l r f + V • r f v (5 l + ■ T 


,3 , 3 

ii + r f 


- (Ii - j) 3 - (Si + j) : 


(76) 


The mass of the body at time t is then given by 


m(t) = m Q - pAV vap , 


(77) 


where m Q is the initial mass of the body. 


VII. PHYSICAL PROPERTIES OF THE BODY MATERIAL 


Physical properties of the glassy material pertinent to the abla- 
tion problem are the following: 

(a) The thermal conductivity, k. (kcal/m°K sec). 

(b) The density, p. (kg/m 3 ). 

(c) The specific heat at constant pressure, Cp. (kcal/kg°K). 

(d) The emissivity constant at the surface, e. (-) 

(e) The molecular weight of the vaporized gas, M va p. (kg/kg mole). 

(f) The heat of vaporization, hy. (kcal/kg). 

(g) The viscosity, p (kg/m sec) represented by the function 


p = Bi exp 


Bp 


T - B- 


t » 


+ B 4 , (kg/m sec) 


where Bi, B 2 , B 3 , and B 4 are constants. 


(78) 
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(h) The vapor pressure P va p (kg/m sec 2 ) represented by the follow- 
ing function recommended by Chapman [3] which accounts for vaporization 
suppression by oxygen: 


vap 


(P /P* ) 

' e vap 


m n 


(kg/m sec 2 ) 


(79) 


where mj, = constant. The equilibrium vapor pressure of the vaporized 
gas P* is given by 


P* 

vap 


= A x exp 



(kg/m sec 2 ) 


(80) 


where Aj., A 2 , A 3 , and A 4 are constants. 

(i) The refractive index, n. (-) 

(j) The absorption coefficient, a^. (1/m). 

(k) The reciprocal radiation mean free path, a. (1/m). 

(l) The effective reflectivity of the surface, R e ff. (-). 
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APPENDIX A 


The Differential Equations of the Glass Layer 


The three differential equations describing the viscous glass layer 
in the vicinity of the stagnation point are well known. The equations 
of continuity, momentum, and energy, where the reference system is 
fixed at the stagnation point, with independent variables X, Y measured 
in the direction shown in Figure 6 with the corresponding velocity com- 
ponents are as follows: 

Continuity 


u Su Sv 

X + Sx BY 


0 . 


(A-l) 


Momentum 


S f Su\ dP oX dW 

37 3yJ = dX ’ dT * 


(A- 2) 


Ener &Z 


pc 


P 


St . S 2 t 
a?* k 3F' 


St 

pC P V SY 


Sf 

Sy‘ 


(A-3) 


Because of the large Reynolds number of the glass-liquid layer, the 
inertia terms have been emitted, since they are negligible in compari- 
son to the shear and pressure gradient terms. The variation of the 
pressure in the Y direction is assumed to be zero which is consistent 
with usual boundary layer equations. 

In the vicinity of the stagnation point the velocity component u 
varies linearly with X; i.e., u = cX. The continuity equation (A-l) 
can therefore be written as 


X Sv 
2 SY ’ 


(A-4) 
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which when differentiated with respect to Y yields 


du _ X b^v 
bY “ “ 2 bY^ * 


(A- 5) 


Integrating the momentum equation from the surface Y = 0 to an arbitrary 
point in the glass-liquid layer yields: 


where 


5u _ dP 
^ Sy " dX 


0X dw 

dt T w 



(A-6) 


(A- 7) 


Substituting equation (A-5) into (A-6) yields 



2_ dP + 2^ dW 

X dX r f dt_ 


Y + 2 — 
X 


(A- 8) 


The Newtonian pressure distribution yields the pressure term in the 
vicinity of the stagnation point as 


I dP 
X dX 



(A- 9) 


Integration of equation (A-8) twice, first from Y 0 to Y and then from 
zero to Y and due to equation (A-9) and the boundary condition bv/bY = 0 
at Y = Y 0 results in the following equation for the ablation velocity v 


Y Y r 


v(Y,t) = v (t) 
w 


o Y 


rj (t) 


P e (t) -P M (t) 




dW(t)l 
(t) dt J 


Y + 2 


T (t) 
W 

X 


n(Y,t) 


dYdY. 


(A-10) 
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The flux of radiative energy represented by the symbol F which 
appears in the energy equation (A-3) is given by (see Kadanoff [9]) 


2r^a j 




T 4 (q) exp -a(Y - tj) dq - / T 4 (tj) exp -a(T) - Y) dq 


+ R « / T 4 (q) exp -QJ(Y + q) dr\ 


where n is the refractive index, CL A is the absorption coefficient, a i 
the Stefan-Bol tzmann constant, a is the reciprocal radiation mean free 
path, and R e ff the effective reflectivity of the surface. 


Figure 6. Space Variables (X,Y) and Corresponding Velocity 
Variables for Surface Fixed Reference System 





APPENDIX B 


The Fortran Program and Its Input Data 


A. Preparation of the Input Data for the Computer Program 

The computer program and the subroutines employed by the pro- 
gram are presented in Fortran IV language in Part B of this Appendix. 
To run i this program, the following input cards must be prepared: 


Card 

No. Columns 

1 1 


2-3 

4 


5-16 

17-32 

33-48 

49-64 

2 1-16 
17-32 
33-48 
49-64 


Integer denoting the method (Method 1 or Method 2) 
to be used for calculating the front face radius 
and mass of the body. Should be 1 for Method 1 
and 2 for Method 2. II format. 

Not used. 

This column is used to give alternative of whether 
to include internal radiation effects (with internal 
radiation) or assume glass is opaque and account 
for only heat flux radiated away from the surface 
(without internal radiation). Should be 1 for 
case without internal radiation or 2 for case with 
internal radiation. 11 format. 

Not used. 

Initial body shape, sphere or hemisphere; should 
be 0 for hemisphere, non-zero for sphere. 

(thickness grid along axis). 

r f (initial front face radius). 
r o 

H q (initial altitude). 

W Q (initial velocity). 

0 Q (initial entry angle in degrees). 

(altitude where slip flow regime begins, 
see Figure 5). 
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Columns 

1-16 

17-32 

33-48 

49-64 

65-80 

1-16 

17-32 

1-16 

17-32 

33-48 

49-64 

65-80 

1-16 

17-32 

33-48 

49-64 


k (thermal conductivity of body material). 

p (density of body material). 

Cp (specific heat of body material). 

M^p (molecular weight of vapor). 

h v (heat of vaporization of body material). 

Qy (vaporization coefficient). 

(3 X (constant in shear stress relation). 

Ai (constant in equilibrium vapor pressure 
function) . 

A 2 (constant in equilibrium vapor pressure 
function) . 

A 3 (constant in equilibrium vapor pressure 
function) . 

A4 (constant in equilibrium vapor pressure 
function) . 

m 2 (constant in vapor pressure relation). 


Bi 

B 2 

b 3 

b 4 


(constant in 
(constant in 
(constant in 
(constant in 


viscosity function), 
viscosity function), 
viscosity function), 
viscosity function). 


Card 

No. 


Columns 


e (emissivity constant of opaque body material), 
a (reciprocal radiation mean free path). 

«A (absorption coefficient), 
n (refractive index). 

R e ff (effective reflectivity of the surface). 

At 3. (first time step). 

Hii! (first maximum time). 

MPi (print frequency for first interval; 14 
format, right adjusted to column 36). 

9 1-16 At 2 (second time step). 

17-32 Tm 2 (second maximum time). 

33-48 Mp 2 (print frequency for second interval; 14 

format, right adjusted to column 36). 

10 1-16 At 3 (third time step). 

17-32 Tm 3 (third maximum time). 

33-48 Mp 3 (print frequency for third interval; 14 

format, right adjusted to column 36). 

Unless specifically stated otherwise, all fields are formated E16.8. 


8 1-16 
17-32 
33-48 


7 1-16 

17-32 
33-48 
49-64 
65-80 
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B. THE FORTRAN PROGRAM AND SUBROUTINES 


C THE ABLATI0N PR0GR AM WITH INTERNAL RAD I ATI 0N OPTION 

C F0R INTERNAL RADI AT 1 ON, SET C0L. 4 0F 1ST INPUT CARD = 2 

C SET C0L. 4 0F 1ST INPUT CARD = 1 IF N0 INTERNAL RAOIATI0N 

C IF INTERNAL RAD 0PTI0N LS USED, AN EXTRA INPUT CARD 

C IS REQUIRED, WHICH GIVES THE RADI ATI 0N CONSTANTS 

C THIS BECOME THE 9TH INPUT CARD. 

ABSFf X ) = ABS( X ) 

DIMENSION RGI100), RG1 ( 100) ,TD( 1001 ,ZI 100) 

REAL MUE,KM,MSTAR,M22,M1,K, KWA I R , LH , MXTI M, MU WA I R , MU,MVAP, 
1M00,MASS,M0,MUINF,M,MUW 

EQUIVALENCE ! RFO ,RF0 ) , ( V00 , V0O, VO0, VOO) , ( TLI MT I 2 ) . TMLMT ID ) , 

1( MOO, MOO, MOO, MOO), ( P00 , POO , POO , POO ) ,(M0,MO) , ( T 00 , T 00 , TOO, TOO ,T I NF ) 

1 , ( E60SEC , E60 SEC ) , 

2 (PS, PE ) , <TS,TE>, (RS,RH0E,RHOE) , I RH0 INF , RHO I NF ) , ( A00 , A0O, AO0 , AOO ) , 
3(T0,TO) , (VI, VW) 

EQUIVALENCE ( MUWAIR, MUA IR ) 

DIMENSION I PRT ( 3 ) 

DIMENSION TMLMT ( 3) 

DIMENSION TL IMT ( 4) 

DIMENSION RA( 100),RA1( 100) ,DFDY( 100) 

EQUIVALENCE ( TL I MT ( 2 ) , TMLMT ( 1 ) ) 

DIMENSION TDT ( 3 ) , C0FHW ( 5) , C0FKI 6) , T ( 100) , VV ( 100) , 

1TL AST ( 100),D2Y! 100) .DTDY(IOO) ,YI 100) .DTDT(IOO) .DTDTI(IOO) . 

2YB(200),TMU( 200 ) , ARG! 200 ) , ARG1 < 200) 

DIMENSION PR( 15 ) 

DIMENSION TABVI ( 3) »TABV0(3) 

3900 FORM AT (49X33 HM AT ERIAL PROPERTIES AND CONSTANTS///) 

3901 FORM AT (10X25 HK (THERMAL CONDUCTIVITY) =23XE16.8,26X20H (KCAL/M/DEG K 

1EL/SEC) ) 

3932 FORMAT ( 1H09X5HDY = E13.6.8H METERS) 

3955 FORMAT!/////) 

3903 F0RMAT(1H09X19HCP(SPECIFIC HEAT) =29XE 1 6 . 8 , 2 6X1 6H ( KCAL/KG/DEG K) ) 

3902 FORMAT! 1H09X 14HR HO (DENSITY) =34XE 16. 8, 26X10H ! KG/M**3 ) ) 

3925 FORMAT ( 1H09X25HE (EMISSIVITY CONSTANT) =23XE16.8) 

3906 FORMAT! 1H 0 9X18HALFAV(VAP COEFF ) = 30 XE 16. 8,26X1 7H ZERO FOR INFINITY) 

3904 FORMAT! 1H09X23HMVAP ( MOL WT OF VAPOR) =25XE 16. 8 , 26X12H ( KG/KG MOLE)) 

3905 FORMAT! 1H09X 26HHV ( HE AT OF VAPORIZATION) =22XE16.8,26X9H!KCAL/KG) ) 

390? FORMAT! 1H09X5HYB = 43X14, 3H DY35X10H ( ME TER S ) ) 

3908 FORMAT! 1H09X 19HVI SCOSI T Y FUNCTION /1H018X27HMU = B1«EXP(62/(T-B3) 
1+ 84)54X1 OH! KG/ M/SEC ) ) 

3909 FORMAT! 1H018X5HB1 = E 16 . 8, / 1H0 1 8X5HB2 = E16. 8 . / 1H0 1 8X5HB3 = E16.8, 
1/1H018X5HB4 = E16.8) 

3910 FORMAT! 1H09X23HVAP OR PRESSURE FUNCTION) 

3911 FORMAT! 1H018X25HPVAP= PE/( ( PE/PVAPS ) **M1 ) ) 

3912 FORMAT! 1H018X34HPVAPS = A1*E XP ( A2/TW**2 +A3/TW + A4 ) 47X1 3H ( KG/M/ SEC* 
1*2) ) 

3913 FORMAT (1H018X5HM1 = E 16 . 8. / 1H0 1 8X5H A 1 - E 1 6. 8 . / 1H0 1 8X5HA2 = E16.8, 

1 / 1 HO 18X5HA3 = E 1 6. 8 ,./ 1H0 18X5HA4 = E16.8) 

3914 FORM AT (1H050X29H INITIAL BODY GEOMETRY METHOD II,///) 

3915 FORMAT! 1H09X16HB0DY IS A SPHERE) 

3916 FORMAT! 1HO9X2OHB0DY IS A HEMISPHERE) 

3917 FORMAT! 1H09X22HRF0 (INITIAL RADIUS) =26XE16. 8 , 26X8H ( METERS ) ) 

39)8 FORMAT ! 1 HO 9X 22 HM0 (INITIAL MASS) = 26XF1 6. 8 . 26X4H ( KG ) ) 

3919 FORMAT (1H052X29HIN I TIAL TRAJECTORY CONDITIONS) 

3920 FORMAT! 1H018X21HH (FLIGHT ALTITUDE) =1 8XE 1 6. 8 , 26X8H ( ME TERS » ) 

3921 FORMAT! 1H018X21HW (FLIGHT VELOCITY) =1 8XE 16. 8 , 26X7H ( M/SEC ) > 
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3922 FORMAT ( 1H018X23HPHI JANGLE OF ATTACK) -16XE1 6. 8.26X9H ( DEGREES ) ) 

3923 FORMAT! 1H062X9HGRJD SIZE.///) 

3924 FORMAT! 1H09X5HDT = E13.6.16H BETWEEN TIME = E14.6.12H AND TIME = E 

114.6.33H WITH A PRINT FREQUENCY 0F EVERY I4.11H TIME STEPS) 

4941 F0RMATJ///1OX35H THIS RUN WILL TERMINATE AT TIME = E16.8.4H SEC) 

6931 F0RM AT ( 1H09X 17 HAL FV 1 VAP C0EF ) -33X8HINF INIT Y ) ; 

C BEGIN CALCUL AT 1 0NS 

3 969 F0 RMAT (£16.8,E16.8,I4) 

3960 FORMAT! 5E16. 8) 

1 CONTINUE 

CALL SCL0CK ( DATE* CT IME > ESEC * E60SEC ) 

MM- 10 _ . ........ 

TLIMTJ 1)=0. 

PI- 3.1415927 . _ _ 

G-1.4 

IYO-6 

T0-300. 

4444 F0RMAT(I1,2XI1, 12X f 3E16.8) _____ 

READ (5,4444) METH0D, I RAO, SPHERE »DY ,RF0 

C JRAD=2 F0R INTERNAL RAOLATI0N 0PTI0N __ 

C IRAQ =1 F0R CASE WITHOUT RADIATION 

IF (IRAQ) 8010.8011,8010 

8010 IF ( I RAD-2 ) 8013,8013,8011 

C IF IRAD N0T -1 0R 2 SET =1 . _ _ . 

8011 IRAD-1 

8013 CONTINUE _ 

IF ( METHOO-l) 2344,2345,2344 

C IF METHOD NOT EQUAL 1. SET IT EQUAL 2 AND CONTINUE 

C METHOD DENOTES METHOD TO CALCULATE BODY GEOMETRY 

2344 METHOD-2 _ 

2345 CONTINUE 

READ (5,3960) HO. WO. PHI ,H 7 0 

READ (5,3960) K, RHO, CP, MVAP ,HV, ALFV, SE 1 

READ (5,3960) A1,A2.A3.A4,M1.B1.B2.83,B_4 

READ (5,3960) E, ALFR, ALFA, RN.REFF 

SIGMA-.1378E-10 _. 

GO TO (8015, 8016), IRAD 

6016 E-0. . .. 

8015 CONTINUE 

D03971 IKA-1.3 

3971 READ ( 5 ,3969) TDT(IKA), TMLMT ( IKA) ,IPRT( IKA) 

IYZ= IYO+1 _ 

CONST-2. *RN**2 *ALFA»S IGMA 

»LEY-L__ 

HMAX -HO +50000. 

R-8314. 

AVM =0. 

MXT I M— TMLMT ( 1} 

C SELECT LARGEST ENTRY IN TIME TABLE FOR END OF JOB TIME 

DO 3961 1 = 2,3 . 

IF (MXTIM-TMLMTI I) ) 3962,3961,3961 

3962 MXTIM- TMLMT ( I ) 

3961 CONTINUE 

TAWX-O. . - 

PH 10-PHI +P 1/ 1 80. 

M0-4./3. *PI»RF0*»3+RH0 

IF (SPHEREI626, 625,626 

625 MO-MO*. 5 — 

626 CONTINUE 

3944 FORMAT (1H09X8H8ETA1 = 40XE16.8) ......... 

3931 FORMAT ( 1H0////I 
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C READ INTERNAL CL0CK, PRINT DATE AND TINE 

3927 FO RMAT IIH1I __ __ _ 

3997 F0RMAT ( 1H1 , I00X6HT IME , A6, /101X6HDATE ,A6 , // ) 

WRITE (6,3997 ) CTIME, DATE 

SEC * E60SEC 

G0 T 0 ( 8040, 8041J, IRAD 

8040 CONTINUE 

8439 FORMAT < 43X47H THE ABLATION PROGRAM WITH0UT INTERNAL RAD I T I 0N//// ) 
WRITE (6,8439) 

G0 TO 8042 

8041 CONTINUE 

4939 FORMAT (43X48H THE ABLATION PROGRAM WITH INTERNAL RADIATION ///) 
WRITE (6,4939) 


8042 CONTINUE 

WRITE (6,3900) 
WRITE (6,3901) K 



WRITE 

(6,3902) 

RH0 



WRITE 

(6,3903) 

CP 



00 T 0 

(8017,8018), IRAD 


8017 

C0NTINUE 




WRITE 

(6,3925) 

E 


8018 

C0NTINUE 




WRITE (6,3904) MVAP 

WRITE (6,3905) HV 

IF ( ALFV ) 6932,6934,6932 


6934 CONTINUE 

WRITE (6,6931) 
GO TO 6933 


6932 

CONTINUE 



WRITE (6,3906) ALFV 


6933 

C0NTINUE 



WRITE < 6 , 3944 ) SE 1 
WR I TE ( 6, 3907 ) I Y0 
WRITE (6,3908) 

- - - - - 



WRITE 

WRITE 

(6.3909) 8 1, B2 , 63 , 84 

(6.3910) 




WRITE 

WRITE 

(6.3911) 

(6.3912) 


> 


WRITE 
G0 T0 

(6,3913) M1,A1,A2,A3,A4 
(8020,8019), IRAD 




C PRINT INTERNAL RADIATION CONSTANTS 

8019 CONTINUE 

8000 FORMATt 9X8HR ( EFF ) =41XE16.8) 

8001 FORMAT! 1H08K3HN =46XE16.8) 

8003 FORMAT ( 1H08X8HALF = 41XE16. 8, 26X10H( 1/METERS ) ) 

8002 FORMAT ( 1H08X6HALFA= 43XE16. 8. 26X10H( 1/METERS ) ) 

8004 FORMAT* 1H08X7HSIGMA =42XE16. 8 ,26X26H( KCAL/ ( M**2 SEC DEG K**4))) 

8005 FORMAT ( 1 HI 51 X28H INTERNAL RADIATION CONSTANTS///) 

WRITE (6,8005) 

WRITE (6,8000) REFF 
WRITE (6,8001) RN 

WRITE (6.8002) ALFA 

WRITE (6,8003) ALFR 
WRITE (6,8004) SIGMA 

8020 CONTINUE 
WRITE (6,3927) 

WRITE (6,3914) METHOD 

IF (SPHERE >3928.3929. 3928 

3929 WRITE (6,3916) 

.GO TO 3930 
3928 WRITE (6,3915) 
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3930 CONTINUE 

WRITE <6,3917) RFO _ 

WRITE (6,3918) M0 

WRITE (6,3955) 

WRITE (6,3919) 

WRITE (6,3920) HO _ __ _ _ 

WRITE (6,3921) W0 

WRITE (6,3922) _P_HI _ _ _ _ __ 

WRITE (6,3955) 

WRITE (6.3923) 

C SET UP PRINT 0F TIME STEPS AND PRINT INTERVAL. 

_ D0 3111 MMM= 1,2 _ 

IF ( TMLHT ( MMM)— TMLNTIMNM+l ) ) 3111,2222,2222 
3111 CONTINUE 
MMM=MMM+1 

2222 CONTINUE 

DO 3933 IKJ= 1, MMM 

3933 WRITE (6,3924) TDT(IKJ), TMLMT ( I KJ-1 ) ,TMLMT( IKJ ) , IPRTI IKJ ) 
WRI TE( 6, 3932 ) DY 
WRITE (6,4941) MXTIM 
HKM = H70 

H70KM «H70 

KTR=1 

IPRINT=IPRT( 1) 

Q2T0L=.O1 

QT0L = .OOl . 

DY2 = DY*DY 

RKCP * K/(RH0» CP) 

RCP=1./(RH0*CP) 

TST= TDT (1) 

DIV-MM 

DO 909 1=2,3. . _. 

IF (TDTM I )-TST)909,909,908 

■_9 09 T$T *TPT U ) 


909 

C0NTINUE 



IF (TST/DY2- 1„/(RKCP*PI) 1911,911,910 


C 

GRID RATIO VIOLATED, DY WILL BE CALCULATED 


910 

DY2=RKCP*PI*TST 



DY-S0RTIDY2) 



WRITE! 6,912) 



912 FORMAT ( 43H0GRID RATIO VIOLATED,. DY WILL BE CALCULATED) 


911 _ CONT INUE _ 

TINE=0. 

JM* 0. 

A6* 397. 92346E 12 

A.7± 6,3-?£fr 

A8*23812.9 
.... _ A9-3.15 

A10 = 73. 

.All * .53733E-4 
A12 = 1926.2 

A13 * .715 

A14= .705 

A15 = ,055 

A16 = .755 

A17 = .025 

A18 = 14.76303E-7 

A19 * 110. 

COE HW (1)=.96954E— 13 
COE HW(2)=-2.30515E-9 
COFHWI 3) = 2. 03371E-5 



C0FHWI4) =*24019 
C0FHVU 5 ) =-5* 3944 
001=1.0099865 

C02=-2.722132E-6 

Oil =-4.993 4E -3 
C 1 2 = . 83986E-6 
021= 19.39623E-4 
022=-.860322E-7 

03 1 =-13.931 88E-5 

032 = .419534E-8 

041 =4.1 16558E-6 

042 =- .97149E- 10 

051 =-4. 41 709E-8 

052 =• 856262 E-l 2 
026 = .26 

068 = .68 

O0FKU) = 1.0229805E-21 
O0FK ( 2 ) = -.5989437E-17 
O0FK(3)=1. 387368 9E -14 
C0FKI4) =- 1*81199 7E— 11 
00 FK ( 5 ) = .2429734E-7 

O0FK ( 6 ) = . 67 16646E-8 ; 

RFT=RF0 

SI-1. 

PV AP=0 • 

YS-0. 

YSW=0. 

IK I =1 

RF V=RF 0 

IF t SPHERE) 3940, 3941,3940 

3941 THICK =RF0 
G0 T 0 3942 

3940 THICK =2.*RF0 

3942 CONTINUE 

C CALCULATE NUMBER 0F STEPS IN Y DIRECTION T0 TAKE CALCULATIONS 

AT = THICK/DY +1.5 
N= AT 

IF { N-100 ) 1 1 y 10 * 10 

10 N= 1 00 

11 CONTINUE 

D0 12 1=1, N 

T ( I ) =T0 
TL AST { I ) =T0 
DTDT ( I ) =0. 

D2Y( I)=0. 

DTOY(I) =0. 

DFOYU 1=0. 

Y ( I ) = DY* FL0AT ( 1-1) 

12 CONTINUE 
F0 = O • 

FY8=0. 

VI-O. 

KK = 1 

DT =TDT ( 1 ) 

DT 2 = DT * • 5 
KP R I NT = 0 
H = H0 

HN = H _ 

W-W0 

U= W0* O0S IPHI0) 

V= W0 *S I N ( PH ( 0 ) 
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UN=U 

VN=V _ _ 

RSEA=6380000. 

GSEA =9.80665 

MASS= M0 

KR* 1 _____ 

GO T0 500 

C EN D OF PREC_0MPUTE SECTION 

1000 CONTINUE 

C BEGAN CALCULATIONS FOR A NEW TINE 

TMN = TIME - TMLMT IKK I 

IF (KK-3 ) 5051,5051,696 _ 

5051 CONTINUE 

_IF ( TMM) 5052, 84,84 _ 

5052 IF ( ABS( TMM)— DT/4. ) 84H3.13 

C CHANGE TIME STEP. ANO PRINT FREQUENCY 

84 KK*KK+1 

IK 1*1 __ 

DT* TDTIKKI 

IPRINT = IPRTI_KK» __ __ _ __ 

0T2 =DT*.5 

13 CONTINUE 

C CALCULATE NUMBER OF STEPS IN Y DIRECTION TO TAKE CALCULATIONS 

AT*(THICK -YSI/DY" +1.5 __ 

N*AT 

IF CN-100) 1111,111 1, 1110 

1110 N=100 

till CONTINUE 

IF (N-IYZ) 1112,1113,1113 

1114 FORMAT C29H0B0DY MEL TED AWAY. THICKNESS *I4.3H DY) _ 

1112 KEY=2 

WRITE (6,11141 N _ 

GO TO 507 

1113 CONTINUE 

KPRINT *KPRINT+1 

DO 33 1=1, N _ _ _ 

Tm=TLAST( I I+DTDT ( I l*DT 

33 CONTINUE 

T(N)= T(N-l) 

C FOURTH ORDER RUNGE KUTTA PROCEEDURE. KR DESIGNATES THE PASS 

C KR*1, INITIAL PASS 

C KR* 2.1ST PASS ON _ANY TIME STEP 

C KR*3, 2ND PASS ON ANY TIME STEP 

_C KR*4, 3RD PASS ON ANY TIME STEP 

C KR*5,4TH PASS ON ANY TIME STEP 

FK1U = DUPT»DT 

FKIV = DVDT »DT 

TIME* TM +DT2 _ _ 

KR*2 

U= UN ♦ FK1U*.5 

V* VN + FKIV*. 5 

GO TO 500 

1001 FK2U* DUDT »DT 

FK2V * DVDT»DT 

U * UN + FK2U*.5 

V * VN ♦ FK2V-.5 

KR* 3 

GO TO 500 

1002 FK3U* DUDT*DT 

FK3V =DVDT*DT 

TIME* TM+DT 
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U= UN + FK3U 
V= VN + FK3V 

KR = 4 

G0 T0 500 

1003 

FK4U- OUDT *DT 
FK4V- DVOT *DT 



U= UN+.1./6.*IFK1U +2 .* I FK2U+FK3U ) +FK41M 


V= VN+1./6.*(FK1V+ 2 , * ( FK2V +FK3V)+FK4V) 


KR-5 

G0 T 0 500 _ _ _ 


1004 

CONTINUE 


C 

END 0F RUNGE KUTTA 


C 

COMPUTE MELTING VALUES 


100 

CONTINUE 
M22-M00**2 
TW -Till 

' - ■ 


S=H/ 1000 • 

IF ( W-2100i* I 11405,11409,11409 


11405 

IF (M00-U) 11406,11409,11409 


11406 

PS ri-1 • +M22* • 2 



RE 1=PST1**2. 5 
PSP-RE1*P$T1 



G0 T 0 11408 


11409 

CONTINUE 



PSP= C6./5.* M00**2)**3.5 

* ( 6. / l 7. *M00**2-1. ) ) **2. 5 


RE1 = <(6.*M22)/IM22+5*>)**3.5 

• I 6. / 1 7,*M22-1 . > * (1.+M22/5.J > 


1**2.5 


11408 

CONTINUE 



PSID-PSP* P00 



TEIDAL- T00* ( 1.+M22/5. ) 



DWDT=IU*DUDT+V*DVDT)/W 
C0N= .5*RH0*DWDT 
C0- C01 +C02*H 

Cl- C11»C12»H : 

C2- C21+ C22»H 

C3- C31 + C32»H _ 

C4- C41+ C 42 *H 

C5^ C51 ♦ C52*H 

PSPS= l C C IC5*M00+C4)*M00K:3I*M00+C2)*M00+C1>*M00+C0 

PS- PSPS»PS1D 

IF ( M00 -35.11401,1402,1402 

1402 KM- .63 _ 

RH0E = 2O.*RH0INF 
TE = 45.* T00 
G0 T 0 1405 

1401 IF 1W-2100.) 1403. 1403, 1404 

C CALCULATE TE, RH0£,KM AND PE F0R N LESS THAN 2100 M/SEC 

1403 TE = TEIOAL 
RH0E- RH0INF *RE1 
KM =1.05 

PE =PS ID 

G0 T 0 1405 

1404 CONTINUE 

C CALCULATE TE, RH0E.KM ANO PE F0R W GREATER THAN 2100 M/SEC 

D0 = 4.016949 — 1 .49287E-2 *S 
01 = -.895475 +4.51127E-3»S 

D2 = 9.28796E-2 -.54521E-3*S 

03 = -4.746323E-3 +3.0S088E-5 »S 

04 = 1 1. 6 1 1955E-5 -7.96241E-7 *S 

05 = -10. 86105E-7 +7.84415E-9 *S 

TE = TEIOAL* ( ( ( ( (M00*O5+D4)*M00+O3I*M00+O2)*M00+O1 )*M00 +00 ) 
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RERE = (<U{M00*(3.569O58E-8 +3. 77445E-9+S ) + C-.832101E-5 - 
16.34127E-7+S) >*M00>+ <5.J)6022E-4 +3.2325E-5*S))*M00 _ _ 

2+ l- 1. 334415 E-2 -6.96885E-4* S))*M00 + (. 1975914+ 7.249941E-3 

3*Sn*M00 + . 269623- 2*47746E-2 • S 

RH0E= RH0 1 NF » REl*RERE 



H0 = .868 — • 1325E-2 
HI = -3.98E-4 *S -. 
H2 = 1.2245E-4 • S 
H3 = -1.0185E-5 *S 
BU = W/100G. - 5. 

*S 

09081 
+.022684 
- 1 i 637 E-3 




KM = (< BU*H3 ♦ H2)*.BU ♦ H1»«BU ♦ H0 


1405 

C0NTINUE 



C 

R0UTINE T 0 C0MPUTE 

V MELTING 


20903 

HU E - (SORT C T E ) * ( 1* 50541 E-7 )/(!• + A1 9/ TEH *GSE A 




MUW = ( SORT I TW ) *<1.50541E-7IJ/(1.+Ai9/TW)*GSEA 

R I NRE=RH0 INF/RH0E 

C PR0BSTEIN S EQUATION FOR THE SHEAR STRESS 

TAHX0 = .727M W/RFT*( l.-RINREJ ) ** 1 • 5 *< 8./3.+RINRE)**.75 
1*< (IHUE*TW)/(MUH*TE) |**.447 *SQRT <2.*RH0E*MUW *TE/TW) 

IF (_P_VAP) 1809,1811,1809 

1811 TAWX-TAWX0 

G0 T 0 1812 

1809 C0NT INUE 

TAHX=TAHX0*CS4 + < 1 . + ( SE 1*S I) / ( PS/PVAP-1 . ) )) 

1812 CONTINUE 

D0 10017 IP=1,100_ _ _ . 

10017 VVtIP$=0. 

Og 1 001 8 IP=li 2 QQ 

ARGC IP)=0. 

THUJ_IP) = 0. _ 

10018 ARG1 ( I P ) =0* 

IDU„H=0 _ . . 

DYB= OY/DIV 



NK»(N-1)*MM+1 

00 105 J=1,NK . 

AZK = J-l 

YBAR = DYB+AZK 

YB t J ) - YBAR 

CALL LATLUHl I ERR . I DUM, I 0NE . N, 1 4, YBAR, V, T , TT ) 

IF (IERR) 101,102,101 

1 01 C0N TINUE _ . 

WRITE <6,7911 

79 1 F0RMAT ( 1X5OHERR0R IN I NTERP0LAT I 0N PR0CEEDURE F0R SHAL L T GR IP! 
CALL DUMP ( YB AR, YBAR , 1 , T ( 1 ) , T< 99 ) , 1, Y<1) , YC 99 ) , 1 ) 

C ERR0R IN INTERP0LATI0N 

C0NTINUE 

102 CALL SUBMU ( TT,B1, B2,B3,B4,MU,TAG I 

IF (TAG) 104, 103, 104 

103 TMUi J ) ■. 0. 

G0 T 0 106 

1 04 THU( J)=MU 

IF ( J-200 ) 105,103,103 

105 C0NTINUE ■ - - 

JJ*NK 

G0 T0 9918 - 

106 IF € J-l » 107,107,9919 

9919 C0NTINUE 

JJ-J-1 

9916 C0NTINUE - - - 

Y0= YB< J J ) 



IF UJ-4) 107 » 109* 109 

107 0 0 108 JP=1 , N __ 

C SET ARRAYS EQUAL ZERO IF PROFILE IS N0T LONG ENOUGH 

TMU( JP)=0. 

108 VVLJP) *0. 

_ JV00 = 0. 

G0 T 0 116 
109 CONTINUE 
20902 CONTINUE 

COFFA = 2.»(2.»(PS-P00)/RFT»»2 + RHO»DWPT/RFT ) 

COFFB = 2.* TAWX 
DO 20904 IL=UJJ 

20904 AR6(IL) = t C0 FFA*YB ( I L ) + COFFB) /TMU(IL) 

GO T 0 10902 

10902 CONTINUE 

CALL SUBV1 C ARfcl ,DVB,ARG,JJ) 

502 CONTINUE 

10903 CONTINUE 
DVY = ARG1 < J J ) 

DO 2712 IK=1,JJ 

2712 ARG(IK) = ARG1 f I K ) -DVY 

CALL SUBVIt ARG1.DYB»ARG» J J ) 

V00 = AR-G1 ( J J ) 

JVK=JJ/MM+1 
D0 2713 JZ*1 1 JVK 
JV = ( JZ“1 ) *MM +1 
TMU(JZ) = TMU(JV) 

2713 VVtJZ>= ARGUJV) 

IF (JV-JJ) 9796,9795,9795 

C SPECIAL CASE WHERE PROFILE TERMINATED ON EVEN GRID STEP 

9795 VVfcJVK )=V00 

9796 CONTINUE 
JZ-JVK 

504 CONTINUE _ 

JX» N-JZ 

IF ( JX ) 1 16, 1 16, 1 15 

115 CONTINUE 

D0 1115 IZ=1,JX . 

JVI = JZ +IZ 

TMil(JVI) =0« 

1115 VV(JJVI)=V00 
116 CONTINUE 

NN »N-1 
DO 34 1=2, NN 

D2 Y ( I ) =IT(I-1)-2**T(I) +T(I + 1))/DY2 

DTDY(I) = <T( I+l>- T( 1-1) )/<2.*DY) 

DTDTI1I) = RKCP* D2 Y (I) 

34 CONTINUE 
DTDY ( N ) =0 • 

D2Y ( N ) = 0 • 

DELT = ABSF { T < 1 ) - TLAST(l))/4. 

IF ( DELT- « 5) 35, 36# 36 

35 DELT = .5 

36 CONTINUE 
TW= T ( 1 ) 

I TER = 1 

GO TO (8021, 8022 ), IRAD 

6022 CONTINUE 

DO 7098 NNN= 1 , N 

IF ( T(NNN)-310. ) 7099,7097,7097 
7097 CONTINUE 
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7098 C0NTINUE 

NNN=N _ _ 

7099 CONTINUE 

IF (NNN-IYZ) 7096 » 7095 , 7095 

7096 IF (NNN-1) 7092.7092,7091 

_7O?2_F0»O. _ _ _ 

FYB=0. 

G0_ T 0 7093 _ 

7091 NNN= I YZ 

7095 C0NTINUE 

D0 701 L=1 , NNN 

^D0 700 J=1,NNN _ _ _ 

700 RG U) = (T(J> )**4 *EXPI-ALFR*CY<L1+YC JM) 

CALL C INTO (NNN,DY t RG,RGl) _ 

RA1(L) = RG1 ( NNN ) 

NP1=NNN-L+1 

D0 703 J=1,NNN 

MNi =J-L+1 _ 

703 RG ( MNI ) = (T(J) )**4 *EXP(-ALFR*(YC Jl - Y ( L ) ) ) 

CALL CINTD ( NP1 , DY, RG , RG11 

ARG ( L ) = RGKNPl) 

D0 704 J = 1.L 

704 RGtJ)=(T(J) )**4 *EXP (-ALFRM Y(L) — Yf J) )) 

CALL CINTD ( L , DY ,RG, RG1 ) _ _ _ „ _____ 

RA ( L ) = RG1 (LI 

701 C0NTINUE __ __ 

D0 706 1=1, NNN 

_y. 06 OFDY ( I ) = C0 N ST *<2,*T t ?) -ALpR » (RAm +AR GU ) +R £F F 

1*RA1(I))) 

F0 = C0NST * ( REFF -U) » RA1 (I I 

FYB = C0NST * (RA(IYZ) - ARG ( I YZ I ♦ REFF * (RAl(IYZ))) 
NNN=NNN+1 

7093 D0 7094 I=NNN,N 

7094 DFOYd 1=0, 

8021 C0NTINUE 

200 HE = W2/8374.88_ + .24*TIJNF_ 

IF (TW) 8091,8091,8093 

8092 F0RMAT ( 1H11X32HT W NEGATIVE 0N I TER ATI 0N NUMBER 14) _ 

8091 WRITE (6,8092) ITER 

KEY =2 

G0 T 0 507 

80 93 C0 NTINUE 

HW=C0FHW( 1) 

D0_37 1 = 1,4 _ _ 

37 HW»TW* HU * C0FHW(I+1) 

HEW=HE-HW 

QRAD= 1* 378E— 11 • E • TW**4 

_ _QRADG = 0. 

IF (2100. -W) 38, 39,39 

38 _ WC» SQRT ( A6 / ( A7+HJ ) 

QC0NT = A8 *SQRT (RH0INF/RFT)* (W/WC)»*A9« f HE-HW) / ( HE-A10 ) 

G g.T fl 44 : 

39 C0NTINUE 

MUWAIR =SQRT (TW) • A18/(1.+A19/TW) _ 

KWAIR=C0FK(1 ) 

D0 40 1=1,5 _ _ 

40 KWAIR = TW*K WAI R+C0FK ( 1+1) 

TWTS= TW/TS 

IF (1.6-TWTS) 41,41,42 

41 FNU = A16 +A17*TWTS 

G0 T 0 43 


42 FNU -A 14 + A 15*T WTS 

43 C0NTINUE : 

QC^T^~1<WAIR/SQRT«MUWAIR) * ( TS-TW) *SQRT (TS*RS/TW) 

1» FNU* A 13»» »4 *SQRT ( 1 « 05* W/ ( 2« »RFT) ) 

44 C0NTINUE 

G0„T0^ (460t460t.45)rKTR _____ 

460 CONTINUE 

J3M0L = A11*RH0INF*W*(W2 - A12MTW-T00I) 

G0 T0 46 

45 QAER0 = QC0NT 

G0 T 0 47 

46 QSLIP = ( QM0L*QC0NT ) /SQRT <QM0L**2 ♦ QC0NT**2> 

IF (H-H70KM) 619, 619,622 

619 G0 T0 (620, 62 1 1_4 5 ) , K TR 

C KTR = 1 WHEN H IS AB0VE H70 KM, KTR=2 DURING TRANSIT 1 0N, 

C KT ft= 3 AFTER TRANSIT I0N 

620 KT R =2 

JTR2 = T I ME_ ♦ 10.*DT 

TR1 = TIME 
_ G_0 T0 622 

621 IF (TR2- TIME) 624,624*623 

624 KTR=3 

G0 T 0 45 

623 QACRJ_ =QSL1P_ + j QC0NT-QSL IP I * ( TI ME-TR1 )/{TR2-TRl) 

G0 T 0 47 

622_ JMNUNUE 

QAER0 =Q SLIP 

47 C0NTINUE 

Tl* l./TW 

PVAP_=A1* _EXP I ( A2*T1 +A3)*T1 +A4) 

PV APS=PVAP 

PVAP=PVAP**M 1 *PS**< 1.-M1F 
IF (PVAP-PS) 1515,1516,1516 

C C0RRECT GUESS F0R TW : 

1516 TWMTW-TLASTU))* .5 +TLAST<1) 

T(l ) = TW ___ _ 

G0 T 0 200 
1515 C0NTINUE 

IF ( ABSF(PVAP) - l.E-20) 48,48,49 

48 SI»1« 

PVAP=0. 

CWFQ=0. 

VI*0. 

G0 T 0 52 

49 CWEQ = 1./I1.+ M/MVAP • (PS/PVAP- 1.)) 

SI =il._-CMEQ ) / ( 1 »~CW£Q*( 1>* C68»< M/MVAP) »»C26) ) 

IF <ABS(SI-i.) -l.E-7) 48,48,9797 
9797 CONTINUE 

IF ( ALFV ) 51,50,51 

50 VW = -I • /RH0 *SI *QAER0/HEW • CWEQ/i U-CWEQ) 

G0 T 0 52 

51 AVM =SQRT <2«»PI »R • MVAP • TW ) ✓ ( ALF V* (( PS-PVAP ) »M 

1+ PV AP*M VAP ) ) 

VI *1 •/ ( 2* *AVM*RH0 ) • 1 1«— CWEQ+AVM* SI*QAER0/HEW - 
1 SORT (<1.-CWEQ+AVM*SI*QAER0/HEW)** 2 +4. *A VM*C WEQ* SI 
2*Q AER0/HE W ) J 

52 CONTINUE 

Q2* QAERO • SI + RH0» VI*HV 

G0 TO (8077, 8088), IRAD 
8077 * Q2*Q2-QRAD 
8088 CONTINUE 
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8111 

Q22=-Q2 


8113 

DTDY C 1 ) =Q22/K 



00 211 IQ=lt N 


211 

DTDT (IQ) =DTQTI ( IQ ) -DTDY( IQ)*(VVI IOH-VW) 



1-1./(RH0*CP)*DFDYC IQ) 



T ( 2 » = .33333333 * (TW-T(4)) + T I3» 


1560 

IF I T( 2) -T0» 1560,1561,1561 
T( 2 ) =T0 


1561 

C0NTINUE 



DTDT ( N ) =0* 



D0 53 1=1,2 
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DTOTCI) = C T ( I ) -TLASTtI))/DT 



DT OY ( 2 ) =1 TC 3) -T( 1 ) )/( 2.*DYI 



DTDYC3) = (T44)— T42))/f 2. *DY) 



D2Y(1)=1./K»( RH0»CP*( DTDT( 1 ) +VI*DTDY( 1 ) ) +OFDYC 1 ) ) 
D2Y(2)=1./K»(RH0«CP»(DTOT(2)+(V1+VVI2> )*DTDY(2) ) 


1+DFDY ( 2 ) ) 


- 

02Y ( 3) = ( T { 4) -2.*T<3) +T(2))/DY2 

DTDT ( 3 1 =RKCP*D2Y( 3 )-DTDY( 3)*(VV(3)+VW) 

. 


1-1./(RH0*CP)*DFDY(3) 



D0 54 1=1, IYZ 


54 

ARG( I ) =DTDT ( D+DTDYII)* (VV(I)+VW» 



CALL CINT (ARG.DY.FI, IYZ) 

Q1 = -K »DTDY (IYZ) + RH0*CP * FI 


556 

Q1-Q1+FYB-F0 

C0NTINUE 


506 

CONTINUE 

IF(,ABS(Q2)-QT0t) 65.65*579 


C 

579 

IF ABS(Q2) IS LESS THAN QT0L, THE I TERATI0N 
C0NTINUE 

IS SKIPPED 


T0L=ABS(Q2*Q2T0L ) 
D=Q2-Q1 



IFIABSF(D) - T0L )65,65,55 


55 

IF (ITER-2) 56,60,60 


56 

E2*D 



TW2=TW 



IF (0) 57*57*58 


57 

TW =;TW-QELT 



G0 T 0 61 


58 

TW= TW+DELT 


59 

G0 T 0 61 


60 

El* E2 


TW1=TW2 
E2 = D 

TW2 =TW 

IF ( ABSF( E2-E11-. 000001 ) 65.62.62 

62 

TW= (TW1*E2 - TH2*E1)/IE2-E1) 
AT S= TW-TW2 



IF (ABS(ATS)-IOO.) 9771,9772,9772 


C 

LIMIT THE CHANGE IN TEMPERATURE T0 100 DEGREES MAX. 

9772 

TW= TW2 + 100.*ABS( ATS ) /ATS 


9771 

C0NTINUE 



IF (ITER-15) 61,61,1661 


1661 

C0NTINUE 


712 

F0RMATI46H1TERATI0N F0R TW FAILED T0 C0NVERS 
WRITE (6,712) 

E.END 0F RUN) 


KEY *2 
G0 T 0 507 


61 

ITER=ITER+1 



T ( 1 ) = T W 



C0NTINUE 



G0 T0 200 
65 CONTINUE 
T(l)=TW 

C ITERATION CONVERGED FOR TW 

V00=V00+VW 
00 650 1=1, M 
650 vv(i»=vvm +VW 

C SET UP TABLES 0F VOO ANO VW T0 INTEGRATE 0VER TIME. 

TABVI(3)=TABVI(2» 

TABVII2>=TABVII I) . 

TABV0(3)=TABV0I 2) 

TABV0(2)=TABV0m 
TABVI(l) = VW 
TABV0I1J = V00 
IF ( IK 1-2 1 752,653,654 

752 YS = YS-DT»TABV0( 1) 

YSW=YSW+DT*TABVI (1) 

G0 T 0 655 

653 YS*YS-OT2*(TABV0(I)+TABV0(2) ) 

YSW=YSW+DT2*<TABVI ( D+TABVI (2)) 

G0 T 0 655 

654 YS»YS-DT/12.«(-TABV0(3)+8.«TABV0<2>+5.»TABV0(1) I 

YSW=YSW+DT/12.*(-TABVI< 3 ) +8 . *TAB VII 2 ) +5 . *TAB VI (1) ) 

655 I K I = IK I i-,1 

IF (ABS( YSI-l.E-6) 66,66,9707 
9707 CONTINUE 

GO T 0 (161,162),METH0O 

161 CONTINUE 

RF0YS = RF0-YS 
IF ( RF0YS ) 6796,6796,16111 
6796 CONTINUE 

IF ( SPHERE ) 6798,67,6798 
6798 CONTINUE 

WRITE (6,6797) 

6797 FORMAT ( 1X25H HALF OF BODY BURNED AWAY ) 

KE Y=2 

GO TO 507 

16111 RFT = RFO ♦ .5 • YS**2/RF0YS 

MASS = MO -RH0»PI*YS*(RF0*(RF0-YS*.5)+YS**2/6.) 

RF V = RFT 

GO TO 163 

162 CONTINUE 

RFT = RF0*tl. + .5* ( l.-EXP ( -10. *YS/RF0 I) 32 * < YS/RFO » **2 ) 
LH = RFT -RFO +YS 

Z1 = ( RF0**2 - RFT»*2 + LH**2 ) / ( 2 .*LH) 

ETA1 = (RF0«»2 -Z 1»»2 ) 

IF IETA1) 10904,6641,6641 
10904 CONTINUE 

ETA1=RF0**2 
GO T 0 6641 

6641 ETA1=SQRT(ETA1I 

FC1 = RFO -*• YSW 

RFV=.5*(Z1+FC1) + ETA1**2/(2.*(Z1+FC1>» 

FJS =RFV - RFO -YSW 

MASS =M0 -RH0*PI*CRF0**2 * ( Z 1+RF0 I- 1 . / 3. * ( Z 1** 3+RF0*»3 ) 
l-RFV*»2»(Zl-*-FCl) +1./3. * ( ( Z 1-F JS > **3 ♦ ( FC 1 + F JS )**3 ) » 

163 CONTINUE 

IF ITHICK-ABSF(YS) >67,67.66 

713 FORMAT (28H1B0DY BURNED AWAY, END OF RUN ) 

67 WRITE 16,713) 

KEY =2 
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G0 T 0 507 

66 C0NTINUE 

IF ( KPRINT— IPRINT ) 508,7705,7705 

7705 CONTINUE 

KEY=1 

5O7C0NT INUE 

C KEY * 1 IF REGULAR PRINT, IF TERMINAL PRINT KEY =2 

PHIT= 180. /PI • ATA N2( V ,U) 

YSN=-YSW 

G0 T0 <8023. 8024). IRAD 

8026 C0NTINUE 

_ 1711 F0RM AT i 7H 1 T I ME * _E 1 6,8 , /8H0H E16.8, 8H WINF _E 1 6 . 8.BH MACH E_ 

116.8, 8H PHI E16.8.8H DUOT E16.8,/8H UINF E16.8.8H VINF El 

__ 26,8, 8H OU/DT E1 6.8 .8H DV/DT E16.8,8H GO _ E16 .8./8H RFT E16 

3.8,8H MASS E16.8,8H YS E16.8.8H YSM E16.8.8H RINF E16.8 

6./ 8H TINF EI6,8t 8H PIN F E16.8.8H AINF E16.8.8H MUINF E16.8, 

58H GINF E16.8./8H RE E16.6,8H RH0E E16.8,8H TE E16.8.8 

6H PE ... E 16.8, 8H HUE .. E16.8,/8H KM _ E1 6.8 .8 H H E16.8.8H 

7 HE E 16. 8, 8H HU E16.8.8H PVAPS E16.8,/8H PVAP E16.8,8H 

8AVM E16.8,8H QAR0 _ E16.8,8H $1 _ . E16.8.8H F(0) E16.8./8H T 

9AW E16.8, 8H INT E16.8.8H F(YB) E16.8,///) 

WRITE 16,17111 T IME.H, M ,M00,PHIT,DUOT«U« V, DUOT «DVDT ,CD, RFT, MASS, 

1YS,YSN,RH0INF,T00,P00,A00,MUINF,GT,RE,RH0E,TE,PE,MUE,KM,M,HE,HU, 
2PVAPS»PVAP» AVM.Q AER0 .SI .F0 _,TAWXjF I.FYB 

1712 F0RMATUH V* E13.6.4H T= E16.6,5H TP= E16.6,6H TP2= E14.6.5H DF* E 

116.8,4H V=_E13.6,7H DTD T= E14.6) _ 

J=1 

QS 9121 1=1, N 

J=J+1 

1713 WR I T E (6,171 21 YUl.Tm,DT DY(II ,02 Y( I ) ,DFDY ( I ) , V V( I ) , OTOT ( I ) 

IF IIYZ-H 9736,9738,9738 

9738 CONTINUE _ 

IF (J-45I9722.9723.9723 

9723 WRITE (6,97131 

J=0 

9722 C0NTINUE 

9721 C0NTINUE __ 

9737 C0NTINUE 

G0 T 0 8025 

8023 C0NTINUE 

8031 F0R MAT 1 7H1TIME S E 16. 8./8H0H _ E16.8 .8H UINF E16.8.8H MACH E_ 

116. 8, 8H PHI E16.8,8H OMDT E16.8,/8H UINF E16.8.8H VINF El 

26 .8.8H DU/DT £16.8^ DV/OT _E16.8_,8H CD_ E16.8./8H RFT E16 

3. 8 , 8H MASS E16.8.8H YS E16.8.8H YSM E16.8.8H RINF E16.8 

4./8H TINF E16.8.8H PINF E16.8.8H AINF E16.8.8H MUINF E16.8. 

58H GINF E16.8./8H RE E16.8,8H RH0E E16.8.8H TE E16.8,8 

6H PE E16.8, 8H MUE E16.8,/8H KM E16.8.8H M E16.8.8H 

7 HE E16. 8, 8H HU E16.8,8H PVAPS E16.8./8H PVAP E16.8,8H 

BAVM E16.8 ,8H QAR0 E16.8.8H SI E16.8,8H QR.ACL E16.8./8H T 

9AW E16. 8, 8H INT E16.8,///) 

WRITE (6.80311 T IME. H. U.M00.PHIT. DUOT. U.V.DUDT.DVDT. CO. RFT. MASS. 

1YS,YSN,RH0INF,T00,P00,A00,MUINF,GT,RE,RH0E,TE,PE,MUE,KM,M,HE,HM, 

_ 2PVAPS.PVAP, AVM,QAER0,SI,QRAO,TAWX,FI : 

8032 F0RMAT ( 4H Y= E13.6.4H T= E14.6,5H TP= E14.6.6H TP2= E14.6.5H MU= E 

. 114, 7, 4H V= E13.6.7H DTDT* E14.6) 

J=1 

D 0 97 11 1 = 1. N 

J=J + 1 

8033 WRITE (6 , 80321_ YJU .T U I .DT DY ( 1 1,02 Yi I) , TMOII I . VV ( I) ,DTDT( I) 

IF (IYZ-II 9726,9728,9728 
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9726 IF (fm-301.1 9727,9728,9728 

9728 CONT INUE _ __ 

IF (J-45) 9714,9712,9712 

9712 WRITE (6,9713) 

J = 0 

9 713 F ORMAT ( 1H1 ) _ 

9714 CONTINUE 

9711 CONTINUE _ _ ._ 

9727 C0NTINUE 

8025 C0NTINUE 

IF (KEY-2) 9781,9782,9782 

9783 FORMAT ( 1H1, 

11H 6HZETA1 E 16* 8 ,6H ETA1 E16.8.6H H E16.8,6H RFV E16.8.6H J 
1_ E 1 6 . 8 ) _ _ 

9782 WRITE (6,9783) Z 1 ,ET A1 , LH, RFV,F JS 

9781 C0NTINUE 

C THIS PARTS COMPUTES ELAPSED TIME 0F CASE AND PRINTS 

G0 T 0 (8121,8122,8121) »KEY _ _ .. .... ....... 

8122 CONTINUE 

CALL SCL0CKI DATE »CT I ME, E.SEC , E60SEC) _ _ ... 

IF (DATE >8123,8121,8123 

8123 CONTINUE 

TOT.T I M=E 60S EC-SEC 

ISEC1 = T0T_TIM _ _ 

ISEC = I SEC 1/60 

ISEC2 = I SEC_ * 60 _ _ _ 

I SEC 60 = I SEC 1 - ISEC 2 

ISEC3 = I SEC / 60 

ISEC4 = ISEC 3*60 

ISEC5=ISEC-ISEC4 

8124 F0RMAK31HOELAPSED TIME ON THIS CASE WAS 14, 9H MINUTES, 14, 5H AND 
114, 12H 60TH SECOND) 

WR ITE ( 6, 8124 ) ISEC3, ISEC5.ISEC60 

8121 CONTINUE ; 

C END OF ELAPSED TIME PART 

GiLJO (508, 1,506), KEY _ _ 

508 CONTINUE 

LMH_ = TIME - MXTIM ... 

IF ( TMM ) 5049,696,696 

5049 IF (A8S(TMM)-DT/4.) 696,697,697 

696 CONTINUE 

KEY= 2 _ _ 

4902 FORMAT ( 5 1HOTRAJECT0RY TERMINATED BECAUSE MAXIMUM TIME REACHED) 
WRITE (6,4902) 

GO TO 507 

697 CONTINUE 

TM* TIME 

UN =U . _. 

VN = V 

HN*H _ . .._ 

DO 666 1=1, N 

666 TLAST ( I ) =T( I ) 

IF ( KPRI NT-I PRINT ) 515,516,516 
516 KPRI NT =0 
515 CONTINUE 

GO T 0 1000 
500 CONTINUE 

C ROUTINE TO COMPUTE THE TRAJECTORY EQUATIONS 

H=HN+V*( TIME-TM ) 

IF (H-HMAX) 4917,4915,4915 
4915 WRITE (6,4914 ) 


56 



c 

CASE IS TERMINATED IF ALTITUDE EXCEEDS HMAX 


c 

THIS IS NECESSARY T0 PROTECT AGAINST A BOUNCING BODY THAT 

c 

DOES NOT RE-ENTER THE EARTHS ATMOSPHERE. 


4914 

F0RMAT ( 1H13OH80UNCI NG B0DY, J0B TERMINATED.) 


KE Y=2 
G0 T0 507 

4917 

CONTINUE 

IF (H) 4900,4900,4907 


4900 

KEY=2 


4901 

FORMAT (51H0TRAJECT0RY TERMINATED BECAUSE OF IMPACT WITH 

EARTH) 

WRITE (6,4901) 
G0 T 0 507 

4907 

CONTINUE 


C 

USE ALTITUDE ROUTINE TO OBTAIN FUNCTIONS OF ALTITUDE 


PRC1 ) =H 
ERft«0 


CALL PRA63(PR,ERR) 

IF (ERR -U) 1107*1108,1107 


1108 

1109 

WR1TEC6, 1109) 

FORMAT ( 45HOERR0R ALTITUDE FUNCTION ROUTINE, END OF JOB. ) 


1107 

CONTINUE 

POO = PR(_2)* 10000. 



T0i = PRO) 
RH0INF = PRC 6) 



MUINF = PRC 7 ) 
A00 = PR 19) 



M * PR ( 10 ) 

IF (M) 1859.1860,1859 


1860 

M* 2 8. 9644 


1859 

CONTINUE 



GT *GSEA/< ( H +RSEA)/RSEA)**2 
W2» U**2 +V»*2 



W * SORT CM2) 

Nit * W/A00 


C 

RE s(RH0INF»W *2.* RFOI/MUINF 
CALCULATE DRAG 



IF (H-HKM) 14,14,15 


15 

IF ( M00-9. ) 16, 16 *17 


17 

CD*2 • 

G0 T 0 30 


16 

CD -1.68 + 2.85/M00 
GO TO 30 


C 

IF SPHERE IS ZERO, TEXT I TE IS HEMISPHERICAL, OTHERWISE IT 

IS SPHERICAL 

14 

IF ( SPHERE >19.26, 19 


19 

IF CH00-2.I21.21.2O 


20 

CD = .9 +M00/SORT (RE) 


21 

GO TO 30 

IF (M00- U26)23_t23,22 ..... ..... ... . 


22 

CO * 1.034 -.O27*H00 

M 10. 30 


23 

IF CM00-. 8124,24, 25 


^24 

CD =.5 



GO TO 30 


23 

£D» .812»M00 -.023 


Z6 

GO TO 30 

IF (MOO-2.) 28, ?8.27 


27 

CD-I .35 + M00/SQRT (RE) 
GO TO 30 


28 

CD* 1.35 


30 

C0NTINHF 


TERM2 =( ( .5*RH0INF*CD*W2 +.75*RH0*VW**2) »PI *RF0**2 > / ( W*MASS ) 
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31 

C 

511 


DUDT = -U*V/ ( H+RSEA ) -U*TERM2 

DVDT = -GT + U**2/< <H+RSEA>) -V*TERM2_ 

C0NTINUE 

END 0F R0UTINE T0 C0MPUTE D.E F0R TRAJECT0RY 
C0NTINUE 

G0 T 0 ( lOOOf lOOlj 1002* 1003* 10041 »KR 
END 
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SUBR0UTINE CINTO <N,DEL *GAM,T) 

DI MENS 1 0N GAM IIOOI.TIIOO) 

IFI.N-1) 81,61,82 

81 T<1)*0. 

R ETURN _ _ _ _ 

82 IFIN-2) 83*83*84 

83 TI1)=0. _ _ 

T 1 2) = ( GAMI 1) +GAM 12) ) *DEL/2. 

RETURN 

84 IFIN-3) 85*85*86 
_85 _T I 1 ) =0 . 

Til ) =( GAMI 1 I +GAM 12 ) ) +DEL/2. 

TI3)=(GAM(1)+4.«GAM(2)+GAM(3) >*DEL/3. 

RETURN 

86 IF1N-4) 87.87*88 

87 TI1>*0. 

Tl 2 I = ( GAM( 1) +GAMI2 ) ) +DEL/2. 

T I 3 ) =< GAMI 1 ) +4.+GAMI 2 ) +GAMI 3 ) )*DEL/3. 
T(4)=3.*(GAM(l)+3.«GAM(2»+3.»GAM(3M-GAM(4))*DEL/8. 

RETURN 

88 Tl 1 > s 0» 

T( 2 1 =OEL» I • 348611 1 11*G AMI 1)+. 897222 222 *G AM I 2 )-• 366666667*GAMC 3 ) + . 1 
147222222*GAMI4)-.026388889*GAMI5>I 
T < 3 > =DEL*I .322222222+GAM ( 1 ) +1. 3777T7778*GAMI 2 ) +.266666667+GAMI 3) + . 

1044444444*GAM«4)-.011U1U1»GAM(5M 

J=N-1 

D0 80 I =.4, J 

80 T ( I )=TII-1)+0EL»I.015277778*GAMII-3I-.102777778*GAMI 1-2) +.63333333 
1 3*GAM (I- 1 ) +. 480 555 556+GAMI 1 1 026388889 * GAMI I +1 ) ) 

T l A+l ) =T I J ) +DEL* I-.026388889+GAM I j-3 »♦. 147222222*GAMI J-2 )-. 3666666 

167«GAMIJ-1)».89 7222 22 2»GAHI J)+. 3486111 11+GAMIJ+l) ) 

RETURN 

END 
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SU8R0UTI NE C INK ARG, H, F I , N ) _ 

DIMENSION ARG(IOO) 

M=N-1 

H24=H/24. 

_D0 1919 1=1, M _ _ 

IF (I-1I 1,1,3 

1 F I = H24* ( 9.*ARG«1)+19.*ARG(2) -5,*ARG(3) «■ ARG ( 4) ) 

G0 T 0 5 

2 F I = F I +:H24»{-ARG(I-1) + 13. > ( ARG( I >+ARG( I + 1 > ) -ARSd + 2)) 

G0 T 0 5 

3 IF ( I-M* 2,4,4 

4 F I =F I ♦. H24*(ARG(I-2)-5.*ARG(I-l) +19.*ARGtI) + 9.*ARG<I+1I) 

5 C0NTINUE _ ..... 

1919 C0NTINUE 

RETURN ; 

END 

C SUBR0UTINE T0 C0MPUTE VISC0SITY MU, TEMP IS TEMPERATURE ,B_1,B2,B3 

C ARE C0FFICIENTS.TAG RETURNS AS ZER0 IF MU REAL LARGE , 0THERW ISE 1 



SUBROUTINE SUBHlHTEH P»Bl»B 2 1 B3f.B4f FUtTAS ) 

A=B2/( TEMP-B3) +B4 

IF <A-50.)2,Ltl 

1 TAG*0. 

G0 T 0 3 _ _ _ 

2 FU*Bl*EXPf A) 

TAG-1# _ _ 

3 C0NTINUE 

RETURN 

C END 0F VISC0SITY SUBR0UTINE 

„END _ _ 



>£ or i -si O w 4 s 


SUBR0UTINE SUBV I(FI»H»Y*N ) 

DINENSI0N FI(2<j0),Y<200> 

Fill) = 0 . 

H24 =H/24. 

IF (N-4 ) 1,2,2 _ 

1 G0 T0 8 

2 00 7 1 = 2, N 

IF ( I — 2 ) 3,3,4 

3 Fill) =H24*19»«Y(I-1) *19. » Y(I)-5.»Ytl+l) »YU+2») 

G0 T0 7 

IF I I-Nl 5,6,5 

Fill) =FI I I — 1 1 +H24« (-YII-2) +13.*(YCI-1) +Y(I») -Y{I + 1>) 
G0 T 0 7 

Fill) =FI(I-l) +H24*(Y(I-3l-5.*YII-2) + 19.*Y(I-l)+9.*Y(m 

C0NTINUE 

RETURN 
00 9 I=>1 * N 
Fill) = 0 . 

RETURN 
END 


62 


C LAGRANGE TABLE L00K UP ANO__HULTIPLE INTE RPOLATION 

C BY T0MMY J. HEINTSCHEL 

C GENERAL ELECTRIC COMPANY . FLIGHT ANALYSIS UNIT 

C FORTRAN IV LANGUAGE 

C •_*_* •*****•*•*•_*•***.*** . . *_* .j 

c 

SUBROUTINE LATLUMUERR* IOUM.M.N. P. ARG, XI ,Y1. A NSI ) 

C I ERR = ERR0R SNITCH 

_C IF ERR0R 0CCURS . I ERR = N0NZER0 

C IF N0 ERR0R 0CCURS , IERR = 0 

C _ _ I DUN * PRESENT TABLE L0CATI0N USED BY SUBROUTINE.. 

C BEFORE ENTERING SUBROUTINE FIRST TIME , PROGRAMMER 

C _ MUST SET IDUM=0 _ __ 

C M NUMBER OF DEPENDENT TABLES 

_C N = NUMBER OF TABULAR POINTS PER TABLE 

C P NUMBER OF POINTS USED FOR EACH INTERPOLATION 

C P EQUAL TO OR LESS THAN 10 ... _ 

C ARG = LOCATION OF ARGUMENT « X) 

C XI = LOCATION OF FIRST VALUE OF INDEPENDENT TABLE 

C Yl = LOCATION OF FIRST VALUE OF DEPENDENT TABLES 

C ANSI * LOCATION AT WHICH THE FIRST ANSWER IS TO BE STORED 

(:••••••••••••••••••••••••••••••• 

DIMENSION Xl(N).tVl(NtM) , ANSL(M) 

INTEGER P 

IF IP .GT. 10) GO TO 25 .... 

IERR=0 

IF (IDUM .EQ. 0) IDUH=1 

DO 10 JMDUM.N 

IF I ARG .GT. XII J)> GO TO 10 ... 

GO TO 20 

10 CONTINUE _ 

IF (ARG .EQ. X1IJI) GO TO 50 

2 5 - I E RR .= 9 9 999 

RETURN 

28 J=J-1 _ 

50 DO 40 Ks.l t M 

40 ANSI I K |=Y1I JjKI . .. 

RETURN 

20 IF I ARG .EQ. XHJI1 GO TO 50 

J=J-1 

IF. (J .LT. 1) GO TO 25 .... 

IF (ARG .LE. XII JM GO TO 20 

J=J+1 _ _ 

I DUM= J 

D0NA1=2.»CX1( JI-ARG) 

D0NA2=2.*(ARG-X1(J-1» » 

_ D0NA3=ABS ( XI ( J ) )+ABS(Xl( J-l) ) _ . _ __ 

IF ( ABS( D0NAI/D0NA3 ) .LE. .0000008) GO TO 50 
IF ( ABS( D0NA2/D0NA3 ) .LE. .0000008) GO TO 28 
IF (P-(P/2)»2 .EQ. 0) GO TO 80 

IF I ABS(DONAl) .LE. ABSID0NA2) ) GO TO 80 

ISTRT= J - (P+.1I/2 
GO TO 70 

80 ISTRT- J - P/2 

70 ISTOPs ISTRT+-P - ! . 

IF ( ISTRT .GE. 1 ) GO TO 90 

LS1RI*-J 

I ST0P= P 

GO TO 100 

90 IF ( ISTOP .LE. N) GO TO 100 
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IStRT* N-P+l 
IS10P- N 




100 00 120 K=1,K 




ANS 1 1 K 1 — 0 . 



- 


PPR0D=1. 

00 130 L = I STRT « IVT0P 
PPR0D=PPR0O*< ARG-Xl(L)) 

PR0D=1 • 

00 140 LL=ISTRT, IST0P 

IF (L *EQ« LL ) G0 T0 140 

PR0D=PR0D* IX 1 ( L ) -XI ( LL ) ) 

140 CONTINUE 

ANSI (K) = ANSI (KH-(Y1(L,K)/PR0D)/I ARG-X1CL II 
130 CONTINUE 

120 ANSI < K MANSI ( K| *PPR0O 

RETURN 

END 


64 


SUBROUTINE PRA63CPR, ERROR) 

DIMENSION PRC15) ,PBC14)_tZI<5),PKC6j5)_ t RH0KC6,3),TKC6, 
1ZBL14) ,TMBC14),LMBC14) .DMBC14) ,TB( 14 ) , HB ( 14 1 
REAL LMB.MB.MWT 


1X400010 
5 ) , VTK C 6» 3) * 1X400020 
1X400030 
1X400040 


C 

C 

C 

C 

C 

C 


ENTER WITH PRC 1) =CURRENT ALTITUDE 
CALCUL ATED TABLE AT RETURN 
PR(1**CURRENT ALtlTUDE Z 

PRC 2 > =PRESSUR£ PRES 

PRC3)*KINETIC TEMPERATURE TEMPK 


1X400050 

1X400060 

1X400070 

1X400080 

1X400090 


c 

PRC 5 ) =N0LECULAR TEMPERATURE 

TEMPM 

1X400110 

c 

PRI 6)=DENSITY 

DENS 

1X400120 

c 

PRI 7 ) = V I SC0SITY 

VISC0S 

1X400130 

c 

PRI 8 J =KI NEHAT IC VISC0SITY 

VISK C=0.0 BEY0NO 90j000.0 METERS) 

1X400140 

c 

PR 1 9 ) =5PEED 0F S0UND 

SPDS0 

1X400150 

_c_ 

■PRUQl^HglECUlAR WEIGHT 

MUT 

1X40.0 16.0 


C PRC 11) *SEA LEVEL PRESSURE PSL 
C PRC 12) ■PRESSURE RAT 1 0 PRAT 

C PRC 13)=DENSITY RAT 10 DR 

C PRC 14) “VISCOSITY RAT 10 VR 

C PR l 15 ) * PRESSURE DIFFERENCE DELP 

ERMR»fl. 


1X400170 

1X400180 

1X400190 

1X400200 

1X400210 

1X400220 


1 Z-PRCl) 

IFCZ. GE.O.. AND. Z.LE. 700000. ) G0 T0 20 
ERR0R=1. 

IFCZ. LT. 0.1 Z=0. . _ 

IFLZ.GT. 700000. ) 2=700000. 

20 N C 1 


1X400230 

1X400240 

1X400250 

1X400260 

1X400270 

IX4002B0 


IF C Z -83004.1 40,30,30 
40 IF C Z -ZI(N) )60, 50,50 
50 N*N+1 

G0 T0 40 _ 

60 Z2 -Z*Z 

U =Z2»Z 


1X400290 

1X400300 

1X400310 

1X400320 

1X400330 

1 X4 . 0 0340. 


Z4 =Z2*Z2 
Z5 =Z2*Z3 
G0 T0 100 

30 IFCZ-90000.J 300j 70, 70 
70 IFCZ-ZBCN)I85,400,80 
80 N=.N+1 


1X400350 

1X400360 

1X400370 

1X400380 

1X400390 

1X400400 


G0 T0 70 

85N=N-1 

G0 T 0 400 

KINETIC TEMPERATURE F0R <O_T0 830041 
100 TENPK=TKCl<tN)+TKC2,N)*Z+TK13,N)*Z2 +TM4,N)*Z3 
16, N ) »Z5 


♦TKC5,N)*Z4 


1X400410 

1X400420 

1X400430 
1X400440 
+TM 1X400450 
1X400460 


IFCZ-28000. ) 120,140,140 
C»* »»+ PRESSURE F0R CO_T0 28000) 

120 PRES= 10. 0000000*EXP ( PK C 1,61 ) +PKC 2,N ) *Z+PKC 3, N ) *Z2+PKC 
_ 1,N)*Z4+PKC6,N)*Z5) 

C****« DENSITY F0R CO T0 28000) 

DENS=C1. 16790729 )»EXP CRH0KC 1 , N) +RH0KC7 . N)*Z+RH0KC 3, N)»Z2 


1X400470 

1X400480 

4,N)*Z3+PKC5 1X400490 
1X400500 
1X400510 
♦RH0KC4. 1X400520 


IN ) *Z3+RH0KC 5 *N) *Z4+RH0K C 6»N ) *Z5) 

IF LZ-10832.1) 130, 130,160 
***** VIRTUAL TEMPERATURE F0R CO T0 12000) 
130 TENPV=C VTK Cl ,N) +VTM 2,N) *Z+VTM 3,N) *Z2 
1 ♦VTKC6,N)*Z5) 

G0 T 0 170 


♦VTK C4, N) *Z3 


1X400530 

1X400540 

1X400550 
♦VIK I 5.NW4IX4Q0560 
1X400570 
IX400SR0 


***** PRESSURE F0R C 28000 T0 83004) 

140 PRES=. 000980665* EXP1PKU,N)+PKC2,N)*Z+PKC3,N)*Z2 
1 PKC5,N)*Z4 ♦PKC 6 , N )*Z5 ) 


1X400590 
♦ PK(4,N) *23+1X400600 
1X400610 
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C *»*»* DENSITY F0R <28000 T0 90000) 1X400620 

150 0ENS =34» 8367 6* I PRES/TEN PK ) IX400630_ 

160 TEMPV=TEMPK ~ 1X400640 

C»«**» VISCOSITY <0 T0 90000) 1X400650 

170 VISC0S=( .000 001 458* SORT (TEMPK*TEMPK*TEMPK) ) / ( TEMPK+110.4 ) 1X400660 

C*»*»» K INEMATIC V I SC0SI TY_£0R (0 T0 90000J _ _ 1X400670 

V I SK=V ISC0S/DENS 1X400680 

£***♦* SPEED 0F S0UNO (O_T0 90000) 1X400690 

SPDS0=2O.O468*5QRT( TEMPV) 1X400700 

MWT»28.9644 ; 1X400710 

TEMPM=TEMPK 1X400720 

C**»** VISC0SITY RATIO <0 T0 700000) 1X400730 

180 VR*VISC0S/. 00001830243 1X400740 

C***»* PRESSURE RAT I0_ (0 T0 700000) _ _ 1X400750 

PRAT=PRES/10. 1701472 " 1X400760 

C»**l** DENSITY RAT 1 0 F0R <0 T0 700000) 1X400770 

DR*0ENS/1. 18354674 1X400780 

C*»*** PRESSURE DIFFERENCE F0R <O_T_0 700000) 1X400790 

0ELP=PSL-PRES " 1X400800 

C»*«CALCULATI0NS COMPLETE , RETURN T 0 MAIN PROGRAM _ 1X400810 

G0 T0 500 ' 1X400820 

300 TEMPK*180.65 1X400830 

C***»* PRESSURE F0R 183004 T0 90000) 1X400840 

PRES s PBASE*E XP( (— 1.3733 0 152 3E 12* < Z -83004. ) )/C 180. 65*16344860. +2 1X400850 
1) *(6344860. +83004.) )) 1X400860 

G_0_T0 150 JX400870 

C***»* M0LECUL AR WEI6HT F0R (90000 T0 700000) 1X400880 

400 MWT.=MB(N)+DMB(N) *( Z -ZB(N)) 1X400890 

C****» M0LECULAR TEMPERATURE F0R (90000 T0 700000) 1X400900 

TEMPM=TMB(N)+LMB(N)*(Z-ZB(N) ) JX4009_10 

C***** KINETIC TEMPERATURE F0R (90000 T0 700000) 1X400920 

TEMPK= ( MWT /2 8.9644) *TEMPM _ _ 1X400930 

PRES=EXP ( AL0G( PB (N ) ) + ( 1 . 37330152 3E12/( LMB( N) » ( 6344860. + Z ) *(63448601 X400940 

1 »+ZB( N ) ) ). ) »AL0G( TMB( N) / ( TMB (N) + ( LMB (N)* ( Z -ZB(N)))))) 1X400950 

OENS= 34. 83676*PRES/TEMPM 1X400960 

VISC 0S = ( .000 0014 58*SQRT (TEMPM»TEMPM»TEM PM) )/ ( TEMP M + 110.4) 1X 400970 

V I SK=0 . 1X400980 

VR» V I SC 0S/. 0000 1830243 1 X400 990 

SPDS0=2O. 0468*SQRT ( TEMPM) ‘ 1X401000 

TEMPV=TSMPK 1X401010 

G0 T 0 180 1X401020 

500 PRt2_L=PRES _ _ _ _ 1X401030 

PR(i3) = TEMPK 1X401040 

PR(i4)=TEMPV _ 1X401050 

PROS )=TEMPM 1X401060 

PR ( 6 ) =DENS 1X401070 

PR17 )=VI SC0S 1X401080 

PR ( 8 )=VI SK _ 1X401090 

PR (19 ) =SPDS0 1X401100 

PR610) =MWT _ 1X401 110 

PR (ill) = PSL 1X401120 

PR (-12) =PRAT 1X401130 

PR(.13)=DR 1X401140 

PRtl4)=VR _ _ _ 1X401150 

PR( 15 ) =DELP 1X401160 

RETURN _ 1X401170 

DATA PSL /10. 1701472/ 1X401180 

- DATA PBASE/6.23101759E-5/ 1X401190 

DATA ( ZI ( I ), I = 1. 51/10832.1, 17853.3.28000. ,49000. , 83004./ 1X401200 

DATA ( (PK( 1, JH 1 = 1, 6), J»l r 5)/1. 6871 582E-2 t -l . 14251 76E-4 . -1.361 2327 1X401210 
X E— 9 , 7. 36241^5E— 14, -1 .0800315E-17.3. 30464 32 E- 22 ,-7.9910777E-2,-8.10 1X401220 
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X46438E-5,-5.5522383E-9*3.1116969E-13,-1.6687827E-17,3.8319351E-22,IX40i230 
X9. 8 4142T 7 E-1 ,-2.69 76 91 7E-4 , 8.5227541E- 9 , -3. 9620263Ej- 1 lx 1. Q 1 46471E-1X401240 
X17 *-1.026431 8E— 22» 1X401250 

Xl» 14118495E1 »-4. 11497477E-4, 1. 33664B55E-8, -3 . 59518975E-13, 1X401260 

X5. 10097254E- 18»-2.89055894E-23» 1X401270 

X9. 99 324461, - 2. 58298 177S-4 , 3. 76 139346E-9 , -4. 208 87236S -14, 1X401280 

X1.60182148E-19.-1.92508927E-25/ 1X401290 

DATA C ( RH0M 1»J),1= 1,6) , J=l,3)/1. 3302U 7E-2, -8. 8502 064 E-5, -4.214301X40130 0 
X56E— 9, 5.9517557E-13, -3.9744789E-17, 7.8771273E-22»1 • 2667122E— 1 , 1X401310 

X-l .3373147E-4. 2 « 066737 1E~9» 2. 33961 09E-1 3.-3. 2562503E-17 »7» 9035209E 1X401320 
X-22,9.2751266E-l,-1.4349679E-4,-2.8271736E-9,4.7480092E-14, 1X401330 

X1 .886 3246E-18.-4.2702411E-23/ ____ 1X401340 

DATA I ITM I, J)kI=l,6), < J=l,5)/2.9667877E2,-6.7731001E-3,8.4619805E-IX401350 

_X7i:rl«_70Q404?E—10t l • 1451 454E- 14 1-2« 489878 8£-l9.» 1X401160 

X2.6892151E2, 4.3075352E-3,-‘8.9159672E-7»-2. 8929791E-11,5.0724856E-1 1X401370 

X 5,-1. 1490372 E-19, IX4Q 138P 

X3.706455TE2,— 3.2858965E-2, 2.0645636E-6, -4. 3283944E— 11 .-5.7507242E-1X401390 

X1 7«8.2924583£-21. i _ __ 1X401400 

X2.044798El,2.07698384E-2,-8.63038789E-7,1.66392417E-ll, 1X401410 

X-9430076185E-17.-4.090 051 08E-2 2. 1 X4014 20 

X-4, 9886595 3E 2, 3. 92 13728 IE-2 ,—4.9518060 IE- 7 ,-3.26219854E-12» 1X401430 

X 9.66650364E-17.-4.78844279E-22/ 1X401440 

DATA MVTKII.J), 1=1,6) *J=1»3 ) /2. 9937265E2.-7 . 717628E— 3 * 9. 48672 02 E— 1X401450 

X 7,-1. 7136592 E-10 .1.1 074297E-1 4.-2 .3294094E-1 9, 1 X401 460 

X2.689215 lE2»4.3075352E-3»-8.9159672E-7, -2. 8929791E-11 »5.0724856E-1 1X401470 

X5 , - 1 .1490372E-19. 1X401 480 

X3.7064557E2,-3.2858965E-2,2.0645636E-6,-4.32B3944E-li,-5.7507242E-IX401490 

X17*8.292-4583E-21/ 1X401500 

DATA IZB(l), 1=1,14)/ 1X401510 

X9. E4.1.E5»1*1E5» 1«2E5»1.5E5,1.6E5 ,1.7E5,1.9E5,2.3E5»3.E5»4.E5«5.E51X401520 
X,6. E5* 7. E5/ 1X401530 

DATA (TH B( I ). 1*1. 14)/ 18 0.65, 21 0. 65, 260. 65 , 360, 65. 960.65, 1X4 01540 

X11I0. 65, 1210.65, 1350. 65, 1550. 65, 1830. 65, 2160. 65 .2420. 65, 2590. 65, 1X401550 

X2700.65/ 1X401560 


DATA (LMB( 11,1=1,141/ 3.E-3.5.E-3, 10.E-3,20.E-3,15.E-3,10.E 

XT . E— 3, 5 . E-3,4.E-3.3.3 E -3.2.6E-3. 1.7E -3, 1. 1E-3.1.1E-3/ 

DATA 1*1,14)/ 28.9644,28.88,28.56,28.07,26.92,26.66,26. 

X25. 85,24. 70, 22.66,19.94^17.94, 16.84,16,17/ 

DATA <DMB< I),l*l,14)/ -0.844E-5 .-3.20E-5 ,-4.9E-5,-3.833E-5, 


-3,1X401570 
1X 401580 


40,1X401590 
1X401600 


X2»-2.60E-5.-2.75E-5,-2. 875E-5,-2.914E-5,-2.72E-5,-2.0E-5,-l .1E- 
X-0 * 67E-5 , -0. 67E- 5/ 

DATA «PBI 1). 1*1. 14)/. 1722 4436 IE-4. .31 5 97 17 1 2E-5.. 774389807E-6. 
X.265977111E-6, .535849383E-7, .391284945E-7».295911117E-7 , 

X. 178715656E-7. .739258171E-8 , . 20057 31 16 E - 8. .4304566.06E-9, 

X.117315480E-9, • 37019896 IE- 10, • 128115330E-10/ 

EJU2 


1X401610 
5, 1X401620 


1X401630 

1 X40 1640 


1X401650 

1X40166 0 


1X401670 

1X401680 


*_ 
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